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Abstract 

Coherent  imaging  systems  offer  unique  benefits  to  system  operators  in  terms 
of  resolving  power,  range  gating,  selective  illumination  and  utility  for  applications 
where  passively  illuminated  targets  have  limited  emissivity  or  reflectivity.  In  contrast 
to  incoherent  imaging  systems,  partially  coherent  illumination  causes  difficulty  during 
image  processing  due  to  high  levels  of  image  speckle  caused  by  constructive  and  de¬ 
structive  interference  effects  unique  to  the  highly  coherent  illumination  source.  Image 
speckle  is  caused  by  the  random  phase  delays  that  occur  due  to  target  roughness  and 
the  turbulent  atmosphere  between  the  remote  target  and  optical  system.  To  combat 
such  effects,  a  number  of  short-exposure  images  are  combined  by  incoherent  averaging 
to  arrive  at  an  image  that  has  greatly  decreased  levels  of  speckle.  Unfortunately,  such 
average  images  suffer  from  decreased  spatial  resolution  due  to  blur  resulting  from 
atmospheric  distortion. 

Effective  image  restoration  may  be  realized  by  inverse  filtering  the  recovered 
average  image  with  an  optical  transfer  function  that  describes  the  overall  optical 
system  and  atmospheric  turbulence.  In  cases  where  it  is  inconvenient  or  impossible 
to  measure  the  parameters  of  this  evolving  function,  blind  deconvolution  algorithms 
may  be  applied  to  estimate  both  the  unknown  remote  scene  reflectance,  as  well  as  the 
unknown  system  transfer  function.  This  research  proposes  a  novel  blind  deconvolution 
algorithm  that  is  based  on  a  maximum  a  posteriori  Bayesian  estimator  constructed 
upon  a  physically-based  statistical  model  for  the  intensity  of  the  partially  coherent 
light  at  the  imaging  detector.  The  estimator  is  initially  constructed  using  a  shift- 
invariant  system  model,  and  is  later  extended  to  the  case  of  a  shift-variant  optical 
system  by  the  addition  of  a  transfer  function  term  that  quantifies  optical  blur  for 
given  field-of-views  and  atmospheric  conditions.  The  estimators  are  evaluated  using 
both  synthetically  generated  imagery,  as  well  as  experimentally  collected  image  data 
from  an  outdoor  optical  range. 


IV 


The  research  is  extended  to  consider  the  effects  of  weighted  frame  averaging  for 
the  individual  short-exposure  frames  collected  by  the  imaging  system.  Atmospheric 
distortion  and  laser  speckle  effects  create  difficult  challenges  for  image  registration 
algorithms.  In  addition,  anisoplanatic  image  warping  can  cause  individual  frames  to 
£t  poorly  to  the  aggregate  frame  ensemble.  A  system  is  devised  where  such  frames  are 
automatically  identified  for  removal  from  the  average  image,  and  the  resulting  frame 
average  is  compared  to  the  unweighted  average.  Results  are  presented  to  support  the 
new  algorithm  using  both  simulated  and  experimentally  collected  data. 
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BLIND  DECONVOLUTION  OF  ANISOPLANATIC  IMAGES 


COLLECTED  BY  A  PARTIALLY  COHERENT  IMAGING  SYSTEM 

I.  Introduction 

The  central  focus  of  this  research  is  to  explore  the  challenging  problem  of  image 
reconstruction  of  coherently  formed  images  viewed  by  an  optical  system  with 
a  £eld-of-view  that  often  exceeds  the  isoplanatic  viewing  angle.  The  purpose  of  this 
chapter  is  to  provide  a  brief  introduction  to  the  held  of  image  reconstruction  by  way 
of  blind  image  deconvolution  of  images  obtained  through  a  turbulent  atmosphere, 
and  to  explain  the  particular  difficulties  encountered  with  systems  that  approach  or 
exceed  the  isoplanatic  angle. 

1.1  Speckle  Imaging  Through  Turbulence 

Researchers  have  shown  significant  interest  towards  the  general  problem  of  ob¬ 
taining  accurate  image  estimates  of  a  remotely  viewed  scene  viewed  with  an  optical 
system  imaging  though  atmospheric  turbulence.  A  signihcant  body  of  turbulent  imag¬ 
ing  research  has  been  generated  by  the  astronomical  community,  e.g.  [49,77].  The 
images  obtained  through  these  optical  systems  are  distorted  by  several  effects.  The 
optical  distortions  introduced  by  the  telescope  components  are  fixed  and  relatively 
easy  to  quantify.  A  dramatically  more  difficult  problem  is  the  distortion  induced 
by  the  random  condition  of  the  atmosphere  between  the  telescope  and  the  distant 
star  or  planet.  Additionally,  there  may  be  distortion  in  the  image  caused  by  vibra¬ 
tion  or  motion  of  the  telescope  during  the  integration  period  over  which  the  image  is 
captured. 

Scientists  and  engineers  often  seek  to  deduce  the  degradation  of  the  imaged 
scene  due  to  the  effects  of  a  turbulent  atmosphere.  In  the  case  of  stellar  imaging, 
the  light  from  extremely  distant  stars  travels  undistorted  through  many  light-years 
of  the  vacuum  of  space  prior  to  reaching  the  Earth’s  atmosphere.  There,  pockets 
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of  turbulent  eddies  of  air  with  varying  indices  of  refraction  introduce  random  phase 
delays  on  different  portions  of  the  approaching  optical  plane  wave.  The  net  effect 
is  a  randomly  distorted  speckle  image  formed  by  the  constructive  and  destructive 
combination  of  the  distorted  phase  fronts.  The  effects  of  this  distortion  vary  dra¬ 
matically  in  relation  to  the  length  of  time  allowed  for  image  capture.  Long  exposure 
times  tend  to  average  the  effects  of  the  phase  and  amplitude  variations  to  produce 
a  blurred  image.  In  this  case,  the  average  Optical  Transfer  Function  (OTF)  can  be 
described  as  a  low-pass  hlter,  with  a  cutoff  frequency  dependent  on  the  severity  of 
the  atmospheric  turbulence.  However,  by  reducing  the  exposure  time  to  a  period 
short  enough  to  essentially  freeze  the  motion  of  the  turbulent  media  through  which 
the  plane  wave  must  pass,  a  dramatically  different  effect  is  noted.  In  such  cases,  the 
phase  and  amplitude  distortions  of  the  entire  optical  path  through  the  atmosphere 
tend  to  produce  what  has  come  to  be  known  as  a  speckle  image.  Figure  1.1  shows  a 
simulated  image  of  a  diffraction-limited  point  source  as  viewed  through  the  vacuum 
of  space  without  the  effects  of  a  turbulent  atmosphere.  Figure  1.2  shows  the  same 
point  viewed  over  the  course  of  a  long  integration  period  of  time  through  turbulent 
atmosphere.  The  result  is  a  symmetrically  broadened  image,  and  the  optical  system 
can  essentially  be  regarded  as  having  a  low-pass  OTF  or  broad  Point  Spread  Function 
(PSF).  In  contrast.  Figure  1.3  shows  a  simulated  image  of  the  same  point  source  as 
viewed  through  identical  turbulence  as  in  Fig.  1.2.  However,  this  image  was  obtained 
over  an  integration  period  short  enough  to  capture  the  instantaneous  structure  of  the 
phase  and  amplitude  distortions  of  the  turbulent  media.  The  image  of  Fig.  1.3  clearly 
contains  higher  spatial  frequency  information  than  that  of  Fig.  1.2.  Also  notable  is 
the  global  shift  of  the  image  intensity,  often  referred  to  as  image  tilt  that  results  from 
relatively  large  linear  phase  distortion  components. 

The  overall  average  system  OTF  may  be  regarded  as  the  composition  of  the 
individual  OTFs  that  arise  from  the  hxed  (possibly  aberrated)  optical  system,  the 
turbulent  atmosphere  for  a  given  exposure  time,  and  the  vibration  or  motion  experi¬ 
enced  by  the  optical  system  over  the  same  exposure  period.  For  the  simple  case  of  a 
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Figure  1.1:  Diffraction  limited  point  source  image. 


Figure  1.2:  Long  exposure  average  point  source  image. 


Figure  1.3:  Short  exposure  instantaneous  point  source  image. 
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Linear  Shift-Invariant  (LSI)  system,  the  ensemble  average  system  OTF,  Tisys  can  be 
expressed  as  the  product  of  the  component  OTFs, 


where  Hopt{u,v)  is  the  non-random  OTF  due  to  the  design  of  the  optical  system, 
Hturb{u,v)  is  the  statistically  averaged  OTF  due  to  the  turbulent  atmosphere  over 
some  hxed  integration  time,  and  Hreg{u,  v)  may  be  thought  of  as  the  OTF  formed  by 
the  combination  of  registration  errors  in  the  image  arising  from  vibration  and  other 
linear  motion  components  not  produced  by  the  atmosphere.  Finally,  u  and  v  are 
variables  in  the  spatial  frequency  domain  of  the  image  space.  It  is  important  to  note 
that  approximately  87%  of  the  distortion  caused  by  atmospheric  turbulence  results 
in  linear  phase  plane  tip  and  tilt,  the  effects  of  which  might  be  indistinguishable  from 
translational  motion  caused  by  sensor  platform  motion  and  vibration. 

1.2  Blind  Deconvolution  for  Image  Reconstruction 

In  stark  contrast  to  conventional  deconvolution,  where  accurate  knowledge  of 
the  system  OTF  and  thus  PSF  exists,  the  problem  of  blind  deconvolution  assumes 
that  the  overall  transfer  function  of  the  system  is  unknown.  If  the  system  is  LSI,  then 
the  image  formation  process  may  be  modeled  as 

d{x,  y)  =  o{x,  y)  ®  h{x,  y)  +  n{x,  y)  (1.2) 

where  o  represents  the  true  remote  scene  to  be  estimated,  h  is  the  PSF  of  the  overall 
system,  n  is  additive  noise,  d  is  the  image  captured  by  the  system,  and  ®  represents 
convolution  in  two  dimensions.  The  variables  x  and  y  represent  spatial  coordinates 
in  the  image  plane.  In  many  imaging  applications,  the  noise  is  accurately  modeled 
as  signal  dependent,  often  distributed  as  a  Poisson  random  variable.  Despite  signal 
dependence,  the  noise  process  may  be  represented  as  an  additive  quantity  to  each 
pixel  of  a  formed  image  [4]. 
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The  image  model  given  in  Eqn.  1.2  may  be  strictly  applied  only  to  individual 
frames  collected  by  the  imaging  system.  Many  reasons  might  exist  where  the  system 
operator  requires  more  than  a  single  frame  to  form  a  useful  image.  For  distant  remote 
scenes,  low  signal  levels  might  require  the  summation  of  several  image  frames  to  in¬ 
crease  the  Signal-to-Noise  Ratio  (SNR).  Additionally,  coherent  helds  passing  through 
a  turbulent  atmosphere  often  suffer  from  an  objectionable  degree  of  speckle  noise  due 
to  the  constructive  and  destructive  summation  of  random  phase  fronts  from  individual 
point  sources  that  comprise  the  remote  scene.  As  demonstrated  in  Fig.  1.2,  the  inco¬ 
herent  summation  of  some  quantity  of  these  speckle  images  results  in  a  less  chaotic 
image,  albeit  with  dramatic  attenuation  in  high  spatial  frequency  detail.  For  these 
compelling  reasons,  some  form  of  image  averaging  is  typically  necessary  to  produce  a 
useable  image  for  the  system  operator. 

This  research  focuses  on  remote  scenes  illuminated  by  coherent  light  sources, 
typihed  by  some  realization  of  a  high  peak  power  pulsed  laser  system.  The  incoherent 
summation  of  many  coherent  frames  results  in  an  optical  system  that  may  be  effec¬ 
tively  modeled  as  a  linear  system,  thus  permitting  the  inherently  linear  deconvolution 
operation  in  later  steps  of  image  processing. 

If  h{x,  y)  is  well  understood  and  parameterized,  the  unknown  image  o{x,  y)  may 
be  estimated  using  established  methods  such  as  Wiener  hltering  [10] ,  inverse  hltering, 
recursive  Kalman  hltering,  least-squares  hltering,  and  constrained  iterative  deconvo¬ 
lution  methods  [41].  However,  in  many  cases  of  interest,  h{x,y)  is  also  unknown, 
leading  to  the  body  of  techniques  generally  referred  to  as  blind  image  deconvolution. 

The  held  of  blind  image  deconvolution  is  well  established  in  the  literature  [4, 
13, 14, 20, 21, 37, 43, 46, 48-51,  56, 64-66,  73,  76].  Other  common  terms  describing  the 
technique  include  blind  image  restoration  and  blind  image  reeovery.  A  detailed  pair 
of  excellent  survey  articles  on  the  topic  describes  the  most  promising  techniques  used 
by  image  processing  researchers  [41,42].  The  underlying  assumption  of  this  body  of 
knowledge  is  that  of  linearity  and  shift-invariance  of  the  overall  optical  system.  The 
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problem  of  linearity  was  briefly  introduced  above  and  is  not  considered  a  significant 
issue  given  a  collection  of  incoherently  averaged  short-exposure  image  frames.  The 
problem  of  shift  invariance  becomes  apparent  when  the  optical  held  of  view  begins  to 
exceed  certain  proportions,  and  is  discussed  in  greater  detail  in  the  following  section. 

1.3  Isoplanatic  vs.  Anisoplanatic  Imaging 

Section  1.2  presupposes  some  very  important  and  limiting  properties  of  the  op¬ 
tical  system  used  to  capture  the  images.  Most  importantly,  an  optical  imaging  system 
may  only  be  modeled  using  the  convolution  operation  described  by  Eqn.  1.2  if  it  can 
be  shown  that  the  optical  system  is  linear  and  spatially  invariant.  In  many  practical 
cases,  this  assumption  of  shift  invariance  is  not  valid,  due,  in  part,  to  atmospheric 
disturbances  when  viewing  points  of  a  remote  scene  that  are  separated  by  sufficient 
angle,  but  also  due  to  optical  construction  in  large  aperture  systems  without  atmo¬ 
spheric  turbulence.  The  former  situation  is  of  primary  concern  in  this  research  effort. 

Tactical  sensors  designed  for  use  in  a  battleheld  environment  are  quite  different 
than  those  used  to  observe  distant  astronomical  objects.  A  typical  astronomical 
system  has  a  fairly  small  held  of  view  (FOV),  hence  the  collected  image  may  be 
modeled  by  the  convolution  of  the  remote  object  with  a  single  PSF  [61].  This  PSF 
is  the  Fourier  transform  of  the  average  system  OTF,  T-Csys{u,v),  for  some  choice  of 
long-term  integration  period.  The  length  of  the  integration  time  period  and  other 
details  of  this  statistically  derived  OTF  will  be  deferred  to  the  following  sections. 

Unfortunately,  tactical  sensors  require  a  much  wider  FOV  than  do  astronomical 
telescopes.  Typical  geometry  constraints  of  tactical  sensors  require  that  the  optical 
paths  arising  from  individual  points  that  comprise  an  extended  remote  scene  pass 
through  distinct  parts  of  the  turbulent  atmosphere.  The  system  can  no  longer  be 
well  characterized  as  a  shift-invariant  optical  system,  since  no  single  OTF  may  be 
used  to  describe  the  transformation  of  every  point  in  the  remote  scene  to  the  image 
plane.  An  optical  system  with  a  FOV  that  admits  optical  paths  through  more  than 
one  atmospheric  condition  is  said  to  exceed  the  isoplanatic  angle. 
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Turbulent  Eddies 


Figure  1.4:  Anisoplanatic  imaging  of  point  sources.  Image  paths  through  the  at¬ 
mosphere  are  different  depending  on  the  relative  scene  point  separation.  The  image 
of  point  A  will  be  formed  through  a  considerably  different  atmosphere  relative  to  that 
of  the  image  of  point  B.  Conversely,  images  of  points  B  and  C  will  be  formed  through 
approximately  the  same  turbulent  atmosphere.  The  angle  created  from  the  optical 
axis  to  points  A  and  B  is  said  to  exceed  the  isoplanatic  angle  for  some  level  of  average 
turbulence,  while  the  angle  between  points  B  and  C  lies  within  the  isoplanatic  angle. 


Figure  1.4  depicts  the  geometry  of  a  system  that  experiences  anisoplanatic  ef¬ 
fects.  Paths  traced  from  a  pair  of  point  sources  separated  by  some  distance  to  the 
telescope  aperture  traverse  regions  of  turbulence  that  possess  different  indices  of  re¬ 
fraction  and  thus  tend  to  delay  the  optical  phase  by  varying  amounts.  The  atmo¬ 
spheric  refractive  index  inhomogeneities  or  turbulent  eddies  [31]  are  assumed  frozen 
according  to  Taylor’s  hypothesis  during  the  gating  period  used  to  capture  the  image. 
The  relative  size  of  these  refractive  eddies  causes  varying  levels  of  phase  delay  corre¬ 
lation  between  the  optical  paths,  thus  indicating  a  particular  statistical  structure  of 
the  atmosphere. 

If  the  distance  between  scene  points  is  small,  the  optical  paths  traced  from 
both  points  are  essentially  identical.  In  this  case,  the  transformation  of  the  remote 
scene  to  an  image  behind  the  optical  aperture  may  be  accurately  described  by  a 
single  OTF,  and  the  system  is  said  to  be  spatially  invariant.  However,  if  the  distance 
between  points  is  increased  beyond  some  limit,  the  optical  paths  from  each  point  to 
the  aperture  are  quite  different.  In  fact,  a  separate  OTF  is  required  to  accurately 
describe  the  imaging  transformation  of  each  point  through  the  optical  system. 
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Consider  an  extended  scene  that  can  be  described  as  a  broad  collection  of  points. 
There  exists  some  combination  of  scene  extent  and  atmospheric  condition,  beyond 
which  the  system  may  no  longer  be  accurately  modeled  with  a  single  OTF.  The 
separation  of  the  extreme  points  of  the  scene  gives  rise  to  an  angular  separation  of 
rays  traveling  to  the  aperture.  The  sufficiently  turbulent  optical  condition  is  said  to 
cause  the  system  to  exceed  this  isoplanatic  angle,  hence  the  system  must  be  considered 
spatially  variant.  A  common  definition  of  the  isoplanatic  angle  is  “the  angle  between 
two  points  at  which  their  mean-squared  wavefront  error  due  to  differences  in  the 
atmospheric  path  is  one  radian  squared”  [33]. 

The  isoplanatic  angle  of  an  arbitrary  optical  system  using  spherical  wave  prop¬ 
agation  is  given  by  [61] 


0o{L) 


(1.3) 


where  is  the  atmospheric  structure  constant  ,  A  is  the  mean  optical  wavelength 
and  L  is  the  atmospheric  path  length. 

As  an  example,  for  a  system  viewing  a  scene  at  10  Kilometers  using  a  mean 
optical  wavelength  of  1.54  microns  through  a  nominal  horizontal-path  daytime  atmo¬ 
sphere  with  structure  constant  of  =  10“^"^,  the  calculated  isoplanatic  angle  is  1.1 
microradians.  The  maximum  extent  of  a  remote  scene  is 


d 


max 


2Ltan  • 


(1.4) 


At  a  range  of  10  Kilometers,  the  maximum  spatially-invariant  extent  of  the 
object  under  consideration  is  only  1.1  centimeters.  Most  target  scenes  of  tactical 
interest  will  have  an  extent  that  exceeds  the  isoplanatic  angle  for  moderately  turbulent 
atmospheric  conditions. 

The  ramihcations  of  exceeding  the  isoplanatic  angle  are  significant.  No  longer 
can  simple  linear  deconvolution  be  applied  to  the  images  obtained  from  a  spatially 
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variant  optical  system.  The  relatively  simple  image  model  of  Eqn.  1.2  does  not  apply 
to  such  a  system.  Instead,  Eqn.  1.2  must  be  modihed  to  include  the  effects  of  a  myriad 
of  distinct  PSF  contributions  to  the  model. 

1.4  Previous  Anisoplanatic  Imaging  Research 

Several  research  teams  have  investigated  the  difficult  problem  of  imaging  through 
anisoplanatic  turbulence.  Although  much  of  the  research  has  been  conducted  with 
emphasis  towards  incoherent  imaging  of  celestial  bodies  and  objects  within  Earth’s 
orbit,  there  has  been  limited  research  intended  to  solve  problems  associated  with 
imaging  extended  scenes  across  nearly  horizontal  slant  paths  through  dense  regions  of 
the  Earth’s  atmosphere.  Of  this  limited  horizontal  path  imaging  research,  only  a  small 
subset  has  been  devoted  to  image  reconstruction  using  partially  coherent  illumination 
of  the  remote  scene. 

Roggermann  [61]  has  effectively  applied  a  block-matching  technique  that  treats 
a  captured  incoherent  infrared  image  as  a  series  of  isoplanatic  patches,  each  of  which 
can  be  accurately  modeled  as  a  portion  of  the  scene  transformed  by  a  particular 
OTF.  His  research  team  recognized  that  the  main  effect  of  a  turbulent  atmosphere  is 
to  cause  a  local  linear  phase  delay  or  tip  and  tilt  to  an  isoplanatic  image,  although 
other  effects  such  as  focus  anisoplanatism  occur  to  a  lesser  extent.  In  the  case  of  an 
image  comprised  of  many  isoplanatic  patches,  each  patch  will  undergo  a  certain  level 
of  random  displacement  due  to  the  linear  component  of  phase  distortion  specific  to 
each  patch.  Since  the  propagation  of  an  image  from  the  aperture  to  the  detector  can 
be  approximated  by  a  scaled  Fourier  transform,  this  linear  phase  distortion  causes 
image  displacement  specific  to  each  patch.  The  motion  of  each  isoplanatic  patch  is 
decorrelated  from  the  motion  of  other  patches  in  the  image  to  some  extent.  Given  a 
series  of  independently  realized  images,  a  parallel  processing  algorithm  is  then  used 
to  estimate  the  linear  shifts  experienced  by  each  patch.  The  shifts  are  effectively 
removed  by  the  block-matching  algorithm,  allowing  better  reconstruction  of  the  final 
image  while  retaining  sufficient  high  spatial  frequency.  However,  such  an  approach 
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suffers  from  the  tremendous  computational  burden  of  computing  the  block-match 
derived  motion  estimate  for  each  of  the  patches.  Additionally,  the  algorithm  must 
have  some  prior  knowledge  of  the  turbulence  strength  in  order  to  decide  on  the  number 
of  patches  to  match. 

A  similar  approach  is  employed  by  several  researchers,  although  the  methods 
used  to  align  the  patches  across  multiple  image  frames  has  been  varied.  Fraser  et 
al.  investigated  the  performance  of  a  clever  hierarchical  implementation  of  subimage 
patch  correlation  registration.  Their  theories  were  experimentally  validated  [75]  with 
a  series  of  well  conducted  modelboard  experiments  using  local  heating  elements  to 
cause  optical  turbulence  effects. 

Several  years  later,  Clyde  [15]  realized  good  reconstruction  results  using  gradi¬ 
ent  subimage  registration  techniques  and  found  improvement  over  correlation-based 
methods  reported  in  [23].  A  fairly  comprehensive  study  was  performed  in  [9]  to  eval¬ 
uate  the  effects  of  the  size  of  the  individiual  isoplanatic  patches  required  to  achieve 
acceptable  images  for  application  to  astronomy  and  surveillance. 

Finally,  Bondeau  [6]  derived  a  Bayesian  estimator  to  reconstruct  images  from 
a  series  of  Gaussian  noise  corrupted  edge  contours  presented  to  a  multi-frame  algo¬ 
rithm,  resulting  in  a  reconstructed  edge-map  of  the  scene  with  increased  high  spatial 
frequency  detail.  Essentially,  the  discrete  contour  vertices  compare  to  the  individual 
isoplanatic  patches  described  in  [9,15,23,75]. 

Perhaps  the  most  signihcant  impediment  of  the  application  of  these  and  similar 
algorithms  to  the  tactical  scenario  is  their  poor  performance  in  low  SNR  conditions. 
Given  photon-limited  individual  raw  frames  that  comprise  an  ensemble,  any  registra¬ 
tion  technique  that  must  operate  on  localized  subsets  of  the  entire  image  suffers  from 
relatively  poor  performance  [23]. 

An  innovative  approach  to  recovering  extended  scenes  in  anisoplanatic  imaging 
conditions  is  offered  by  Thelen  [70],  who  uses  phase  diverse  speckle  images  to  jointly 
estimate  the  image  and  parameters  of  several  discrete  phase  screens  used  to  model 
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atmospheric  turbulence.  The  incoming  light  is  split  to  create  the  conventional  image 
together  with  an  image  that  has  a  small  but  known  amount  of  defocus.  Although  the 
estimator  is  not  described,  the  author  claims  to  construct  a  system  that  delivers  a 
maximum  likelihood  estimate  for  the  image,  as  well  as  the  Zernike  basis  coefficients  of 
a  small  number  of  phase  screens  placed  at  various  locations  between  the  remote  scene 
and  aperture.  In  addition  to  estimating  the  original  image,  the  algorithm  allows 
estimation  of  the  component  phase  screens  that  model  the  degrading  atmosphere. 
Such  detailed  atmospheric  information  is  of  interest  to  implementors  of  adaptive  optic 
(AO)  systems,  as  well  as  those  who  seek  accurate  estimate  parameters  describing  the 
structure  of  the  turbulence.  An  interesting  result  of  their  research  is  their  conclusion 
that  phase  screens  more  proximate  to  the  aperture  were  better  estimated  than  those 
closer  to  the  remote  scene. 

A  multiframe  processing  algorithm  is  described  in  [19]  that  has  been  shown 
to  effectively  mitigate  image  degradation  from  coherent  speckle  and  anisoplanatic 
viewing  conditions  by  iteratively  processing  subimage  regions  of  a  remote  scene.  It 
appears  that  the  independent  processing  of  multiple  subimages  by  the  modihed  Ayers- 
Dainty  blind  deconvolution  algorithm  [2]  admits  improvement  for  images  formed  by 
a  spatially  variant  imaging  process. 

A  tributary  of  related  research  is  dedicated  to  the  demonstration  of  the  existence 
of  super-resolution  effects  obtained  by  an  optical  system  that  images  scenes  that 
exceed  the  isoplanatic  angle.  Charnotskii  [12]  postulated  in  1989  the  possibility  of 
achieving  optical  resolution  beyond  the  diffraction  limit  of  a  telescope  by  exploiting 
the  frequency  shifting  components  of  the  turbulent  optical  path  between  the  scene 
and  aperture.  He  then  presented  a  detailed  experimental  procedure  to  observe  this 
effect.  Further  analysis  was  conducted  several  years  later  by  Fried  [24]. 

Gerwe  [28]  devised  an  iterative  algorithm  to  reconstruct  a  remote  extended  scene 
using  a  series  of  short-exposure  images,  and  demonstrated  that  Fourier  components 
above  and  below  the  diffraction  limit  were  enhanced  by  the  technique  [29] .  Addition- 
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ally,  studies  of  the  required  photon  noise  level  were  conducted.  He  later  applied  these 
techniques  to  the  deconvolution  of  under-sampled  images  of  wide  FOV  low-Earth 
orbiting  satellites. 

Most  recently,  Lambert  provides  limited  simulation  data  to  support  the  ap¬ 
plicability  of  this  technique  under  high  signal-to-noise  (SNR)  conditions  [44,45]. 
Horizontal-path  simulations  and  real-world  imagery  are  used  to  support  the  super¬ 
resolution  hypothesis. 

Finally,  it  is  worth  noting  that  Sheppard  has  developed  a  multi-frame  recon¬ 
struction  algorithm  that  apparently  achieves  resolution  beyond  the  diffraction-limited 
cutoff  for  isoplanatic  images,  and  reports  simulation  data  [67,68]  to  support  his  claims. 
Despite  the  postulated  improvements  available  through  super-resolution  techniques, 
the  extremely  high  signal  levels  required  to  achieve  acceptable  imagery  restrict  this 
approach  to  a  fairly  limited  subset  of  the  data  collected  by  tactically  employed  laser- 
vision  systems. 

The  approach  taken  over  the  course  of  the  following  chapters  departs  from  the 
established  body  of  literature  in  several  aspects.  Image  reconstruction  techniques 
that  rely  on  subimage  alignment  suffer  three  major  practical  limitations.  Perhaps 
the  most  fundamental  limitation  is  the  high  SNR  levels  required  to  estimate  the  spa¬ 
tial  displacement  of  each  subimage.  While  correlation  and  block  match  alignment 
methods  have  been  shown  to  work  well  on  large  images,  accurate  alignment  of  small 
subimages  is  only  practical  when  the  imagery  is  relatively  noise-free.  The  choice  of 
the  number  of  subimage  regions  is  also  quite  difficult  and  must  be  based  on  some 
assumption  regarding  the  current  structure  of  the  atmosphere  and  system  FOV.  A 
more  turbulent  atmosphere  would  require  processing  of  many  more  subimage  regions 
than  images  produced  during  relatively  calm  viewing  conditions.  Finally,  the  com¬ 
putational  burden  required  to  align  a  large  number  of  subregions  is  often  in  excess 
of  that  available  on  limited  operational  platforms  in  near  real-time,  especially  under 
turbulent  viewing  conditions. 
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Although  the  modified  Ayers-Dainty  blind  deconvolution  algorithm  provided 
in  [19]  seems  promising,  it  is  not  clear  that  the  convolution  operator  is  appropriate 
for  operation  on  individual  frame  subimages,  since  the  detected  intensity  of  each 
image  is  not  linear  for  coherent  illumination  sources.  In  the  cases  studied  by  Dayton 
[19],  it  appears  that  the  modified  Ayers-Dainty  blind  deconvolution  does  increase 
resolution  for  imagery  presented.  However,  the  application  of  linear  convolution  to 
the  image  restoration  process  for  coherent  imagery  is  not  mathematically  justified  [30]. 
The  approach  of  the  research  in  the  following  chapters  requires  frame  averaging  of 
some  number  of  frames  to  produce  an  average  image.  This  averaged  image  may  be 
accurately  considered  the  result  of  linear  processing  throngh  the  optical  system,  since 
the  incoherently  averaged  image  intensity  at  the  detector  follows  a  linear  relationship 
with  the  intensity  reflectance  of  the  remote  scene. 

The  image  reconstruction  approach  developed  dnring  this  research  is  novel  in 
several  important  ways.  The  estimator  is  developed  nsing  Bayesian  techniqnes  based 
on  the  nnderlying  statistical  model  of  partially  coherent  illnmination.  Althongh  con¬ 
siderable  literature  is  devoted  to  reconstruction  of  incoherent  imagery,  the  approach 
presented  in  this  work  concentrates  on  the  formnlation  of  reconstrnction  algorithms 
specific  to  partially  coherent  illnmination.  The  initial  estimator  is  extended  to  the 
case  where  the  system  FOV  becomes  so  wide  as  to  admit  spatially  variant  effects  in 
the  detected  image.  Rather  than  partition  the  image  into  anisoplanatic  patch  snb- 
regions,  a  transfer  fnnction  is  developed  to  model  the  blur  of  the  entire  image.  This 
approach  is  more  applicable  to  the  imaging  conditions  prevalent  for  tactical  observa¬ 
tion  of  remote  targets  nsing  a  laser  vision  system,  dne  mainly  to  the  low  expected 
signal  levels,  bnt  also  dne  to  the  limited  on-board  processing  capabilities  of  the  car¬ 
riage  platform.  Finally,  a  great  deal  of  emphasis  is  placed  on  the  accnrate  recovery  of 
the  seeing  condition  nnder  which  the  imagery  was  collected.  Snch  an  estimate  may 
be  usefnl  when  developing  imaging  systems  used  for  atmospheric  measurement  where 
scintillometry  techniqnes  become  impractical.  Althongh  the  main  goal  of  the  image 
reconstrnction  process  is  the  formulation  of  usefnl  imagery  to  the  system  operator. 
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accurate  knowledge  of  the  prevalent  seeing  condition  under  anisoplanatic  conditions 
may  be  very  useful  in  some  applications. 

1 . 5  Document  Organization 

The  research  described  within  this  document  is  organized  as  follows.  Chapter  II 
provides  the  necessary  background  and  mathematical  underpinnings  required  to  pose 
the  problem  and  understand  the  fundamental  models  used  to  describe  the  statistics 
of  the  physical  processes  that  occur  in  a  partially  coherent  imaging  system.  That 
chapter  covers  the  models  used  to  describe  the  propagation  of  partially  coherent 
light  through  a  turbulent  atmosphere,  as  well  as  models  to  describe  the  composition 
and  simulation  of  that  medium.  A  model  describing  the  statistics  of  the  partially 
coherent  illumination  source  is  offered  and  explored  in  the  context  of  propagation 
through  a  turbulent  medium.  A  maximum  likelihood  estimator  is  derived  to  establish 
the  free  parameter  of  the  illumination  detection  model,  which  is  used  in  subsequent 
chapters  during  the  application  of  a  joint  estimator  for  the  remotely  imaged  scene 
and  atmospheric  seeing  conditions. 

Chapters  III  and  IV  describe  and  rehne  joint  estimators  based  on  Baysian  esti¬ 
mation  techniques  that  seek  a  useful  solution  to  the  blind  estimation  problem  of  joint 
estimation  of  the  remotely  imaged  scene  together  with  the  seeing  conditions  under 
which  the  data  were  collected.  The  derivations  begin  in  Chap.  Ill  with  a  simple,  shift- 
invariant  model  for  the  imaging  system,  and  are  modihed  in  Chap.  IV  to  include  the 
deleterious  effects  of  imaging  through  anisoplanatic  viewing  conditions.  Simulated 
and  experimentally  gathered  data  are  offered  to  support  the  operational  utility  of  the 
blind  estimation  routines. 

Chapter  V  moves  away  from  the  topic  of  blind  deconvolution  in  order  to  better 
address  the  issue  of  multi-frame  averaging  in  the  context  of  partially  coherent  scene 
illumination.  A  probabilistic  model  is  used  for  the  detected  images  that  comprise 
an  ensemble,  and  this  model  is  extended  to  form  a  likelihood  metric  that  describes 
the  admissibility  of  particular  image  frames  into  the  aggregate  ensemble.  The  chapter 
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begins  by  examining  an  iterative  method  that  assigns  weights  to  each  image  within  the 
ensemble  and  it  is  fonnd  that  the  resulting  weighted  average  image  contains  enhanced 
high  spatial  frequency  content.  The  research  is  continued  in  the  later  half  of  the 
chapter  to  explore  the  feasibility  of  binary  frame  weighting,  whereby  selected  frames 
are  discarded  from  the  ensemble  in  order  to  achieve  similar  increases  in  spatial  detail. 
Simulated  and  experimental  results  are  offered  to  reinforce  the  utility  of  the  binary 
frame  weighting  algorithm. 

The  research  is  concluded  in  Chap.  VI  with  some  remarks  that  demonstrate  the 
applicability  of  the  research  to  several  areas,  as  well  as  recommendations  for  further 
work  in  the  held. 
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II.  Coherent  Imaging  through  Turbulence 

Illumination  of  a  remote  scene  by  a  partially  coherent  light  source  introduces  in¬ 
teresting  possibilities  and  difficult  challenges  to  the  image  reconstruction  process. 
The  body  of  literature  surveyed  in  Chap.  I  is  exclusively  devoted  to  scenes  illuminated 
by  narrowband  incoherent  light.  The  research  documented  in  the  following  chapters 
is  focused  on  the  study  of  methods  to  reconstruct  scenes  illuminated  by  light  that 
is  highly  coherent,  as  might  be  produced  by  a  moderately  stable,  high-power  laser. 
The  Air  Force  and  DoD  components  maintain  great  interest  in  such  systems.  Bene¬ 
fits  include  higher  theoretical  resolution  due  to  the  shorter  wavelengths  compared  to 
forward  looking  infra-red  (FLIR)  systems,  non-reliance  on  ambient  light  conditions 
and  thermal  contrast  ratios,  and  long  range  imaging  due  to  higher  returned  photon- 
count  at  the  imaging  device.  The  high  available  power  levels  from  modern  tactical 
targeting/illuminator  laser  systems,  combined  with  rapid  advances  in  image  collection 
technology  have  made  the  long  range  capability  an  exciting  and  physically  realizable 
feature  of  this  technology. 

This  chapter  presents  background  theory  necessary  to  pose  the  general  problem 
and  conduct  research  toward  restoration  of  an  image  of  a  remote  scene  illuminated 
with  a  coherent  light  source.  Much  of  this  material  is  derived  from  the  established 
literature,  although  several  sections  represent  original  contribution  to  the  held  and 
are  noted  as  such.  Prior  to  the  development  of  a  maximum  a  posteriori  (MAP) 
estimator  for  recovery  of  the  remotely  imaged  scene  in  Chap.  Ill,  several  crucial 
questions  must  be  answered  concerning  the  validity  of  the  models  and  underlying 
assumptions  used  to  construct  such  an  estimator.  The  central  tenants  of  the  sections 
within  this  chapter  are  tied  to  the  fundamental  problem  of  reconstructing  images 
formed  using  coherent  illumination  methods  that  have  passed  through  vast  distances 
of  turbulent  atmosphere.  Strong  emphasis  is  placed  on  the  underlying  statistics  of  the 
physical  imaging  models,  as  well  as  the  random  processes  that  govern  the  turbulence 
between  the  target  and  the  laser  vision  system.  With  a  complete  understanding  of 
the  expected  effects  of  the  atmosphere  on  the  propagated  coherent  target  scene,  the 
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reader  is  provided  the  tools  necessary  to  construct  a  joint  estimator  that  recovers 
MAP  estimates  of  the  remotely  imaged  scene  together  with  the  average  atmospheric 
seeing  conditions  at  the  time  of  image  collection. 

The  following  sections  outline  the  theory  necessary  to  construct  physically  and 
statistically  accurate  models  to  represent  the  propagation  of  coherent  illumination 
through  the  realistic  atmospheric  conditions  typically  encountered  by  tactical  ap¬ 
plications  of  a  coherent  vision  system.  A  model  that  describes  partially  coherent 
illumination  of,  and  reflection  from  a  remote  scene  is  presented  in  Sec.  2.1,  followed 
by  a  brief  treatment  of  the  spatial  sampling  issues  that  arise  during  digital  simulation 
of  such  a  system  in  Sec.  2.2.  The  random  statistics  of  a  turbulent  atmosphere  are 
analyzed  in  the  context  of  creating  accurate  digital  representations  of  turbulence  for 
simulation  of  realistic  remote  images  in  Sec.  2.3.  The  salient  details  of  the  experimen¬ 
tal  imagery  system  used  to  collect  long-range  remote  imagery  data  is  covered  briefly 
in  Sec.  2.4.  The  degree  of  coherence  of  the  optical  illumination  system  used  to  collect 
the  experimental  data  was  not  well  established  at  the  time  of  data  collection.  Because 
of  this,  a  maximum  likelihood  estimator  for  the  speckle  parameter  of  scenes  imaged 
using  a  partially  coherent  system  is  developed  in  Sec.  2.5.  To  better  understand  the 
statistics  of  the  detected  intensity  arriving  at  the  detector  of  the  imaging  camera, 
a  brief  analysis  of  the  statistical  transformations  that  model  the  turbulent  coherent 
imaging  process  are  presented  in  Sec.  2.6.  Some  rather  important  image  intensity 
scaling  and  quantization  effects  on  the  modeled  data  are  discussed  in  Sec.  2.7.  The 
effects  of  registration  and  frame  averaging  are  used  to  justify  a  model  for  the  optical 
transfer  function  imposed  by  the  turbulent  atmosphere  in  Sec.  2.8.  The  chapter  con¬ 
cludes  with  a  discussion  of  the  method  of  knife-edge  OTF  estimation  for  application 
to  seeing  condition  estimation  from  a  series  of  experimentally  collected  laser  radar  im¬ 
ages.  Such  estimates  will  be  used  to  establish  atmospheric  truth  from  experimentally 
collected  coherent  imagery. 
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Figure  2.1:  Simplified  sketch  of  the  optical  path  from  laser  transmitter  to  optical 

detector.  The  remote  scene  is  illuminated  by  a  variable  divergence  laser  illuminator  to 
provide  flood  illumination  of  the  target  scene.  Reflected  light  is  propagated  a  distance 
oi  L  =  zt  to  the  optical  aperture,  and  propagated  again  a  distance  corresponding  to 
the  focal  length  /  to  the  detector  plane. 

2.1  Coherent  Imaging  Model 

Prior  to  exploration  of  the  effects  of  atmospheric  turbulence,  it  is  necessary  to 
construct  a  physical  model  of  the  deterministic  propagation  of  coherent  light  through 
the  atmosphere  between  the  target  and  laser  vision  system.  Figure  2.1  depicts  a 
simplihed  sketch  of  the  imaging  paths.  A  model  used  to  describe  the  formation  of  the 
image  assumes  the  target  is  illuminated  by  a  planar  held 


(2.1) 


with  units  of  volts  per  meter  in  the  plane  of  the  target  a  distance  Zt  meters  from 
the  laser  imaging  system.  The  amplitude  of  the  beam  is  described  by  the  function 
Ab{xn,ym)  and  the  phase  is  Vm)  during  the  coherence  time  r^.  The  variables 

Xn  and  ym  represent  coordinates  in  the  plane  of  the  remote  object  to  be  imaged. 

Although  it  is  tempting  to  make  a  hrm  distinction  between  coherent  and  in¬ 
coherent  illumination,  the  terms  are  actually  extremes  in  a  continuum.  In  practice, 
one  may  obtain  neither  perfectly  coherent  nor  perfectly  incoherent  light.  Rather,  the 
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illumination  must  be  quantified  using  some  measure  of  coherence.  A  long  coherence 
time  indicates  a  highly  coherent  optical  source,  with  the  phase  of  the  optical  waves 
marching  in  lock-step  over  moderate  distances.  In  contrast,  a  source  with  a  short 
coherence  time  will  suffer  from  some  decorrelation  of  the  phase  front  over  relatively 
small  distances.  Another,  perhaps  more  intuitive  way  to  visualize  coherence  time  is 
by  translation  to  the  frequency  domain.  A  source  with  infinitely  long  coherence  time 
will  possess  a  spectrum  that  appears  as  a  Dirac  delta  function,  while  a  short  coherence 
time  source  will  appear  as  a  central  frequency  component,  corrupted  by  noise  side¬ 
bands.  A  narrowband  filtered  incoherent  source  will  have  fiat  power  spectral  density 
over  some  finite  bandwidth  and  will  have  very  short  coherence  time. 

The  incoherent  illumination  treated  in  the  literature  of  Chap.  I  was  actually 
narrowband  filtered  incoherent  light.  Clearly,  such  illumination  has  spatial  and  tem¬ 
poral  correlation,  however,  the  correlation  is  very  limited  due  to  the  relatively  high 
bandwidth  of  the  light.  In  contrast,  the  coherent  illumination  considered  within  this 
research  effort  is  sufficiently  narrowband  that  it  becomes  convenient  to  use  coherence 
time  r  to  describe  its  behavior.  At  some  time  longer  than  the  coherence  time,  the 
phase  relationship  of  the  illumination  is  expected  to  depart  from  that  of  the  reference 
sinusoidal  center  carrier  frequency.  The  coherence  time  of  a  laser  illuminator  may 
be  compared  to  the  integration  time  period  of  the  imaging  detector  used  to  collect 
photons  of  the  illumination  that  reflect  from  a  scene-of-interest.  Detectors  with  rela¬ 
tively  long  integration  periods  or  gate  times  will  collect  photons  over  many  coherence 
periods.  The  significance  of  this  phenomenon  will  be  explored  further  in  this  section. 
Note  that  illuminator  coherence  over  the  duration  of  the  round-trip  travel  time  is  not 
required.  In  fact,  such  long  coherence  times  are  difficult  to  achieve  with  operational 
laser  illuminators. 

The  field  reflected  from  the  target,  u^{xn,ymi  Zt)  with  units  of  volts,  can  be 
computed  by  multiplying  the  field  transmitted  to  the  target  as  described  by  Eqn.  2.1 
by  the  reflectance  of  the  target  r{xn,ym),  times  the  sample  size  employed  in  a  digital 
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representation  A^. 


Vm,  Zt)  =  ^sUl{Xn,  Vm,  Zt)r{Xn,  (2.2) 

where  A  is  the  mean  wavelength  of  the  illuminator  source,  and  9{xn,ym)  is  a  func¬ 
tion  describing  target  surface  roughness  in  meters.  The  held  rehected  by  the  scene 
u^{xn,ym-,  Zt)  is  propagated  back  to  the  receiver  aperture  and  may  be  modeled  using 
a  modihed  Rayleigh-Sommerfeld  integral  designed  to  propagate  radiation  from  point 
sources  [30] 


(^n2  1  ym.2) 


{Xn,yni,Zt)e _ ^ _  J0a(a:„,y™,x„2.y„2) 

^27r  [ztY  +  (Xn  -  Xn2V  +  (Vm  “  ym2)2 


(2.3) 


where  u^{xn2,  ym2)  describes  the  held  at  the  optical  aperture,  and  (paixn,  ym,  Xn2,  ym2) 
is  the  phase  delay  function  caused  by  the  atmosphere,  described  further  in  Sec.  2.3. 
The  held  at  the  detector  plane  of  the  imaging  system  u^{xn3,ym3)  can  be  computed 
using  one  additional  Rayleigh-Sommerfeld  propagation  integral 


Ud{Xn3,ym3)  = 


N2  M2 

E  E 

n2=l  m2=l 


■u^(a;„2,2/m2,0)  e 


j^TT-y/  (/)^  +  (a^n3~^n2)^  +  (yTn.3~^m2)^ 


(/)2  -1-  {Xn3  -  Xn2Y  +  (ymS  “  J/m2)^ 


p3<l>L(Xn'2,Vm2) 


(2.4) 


where  /  is  the  focal  length  of  the  system,  and  is  the  sampling  lattice  spacing  of  the 
optical  aperture.  The  held  u^{xn3,ym3)  is  the  held  in  the  detector  plane  in  units  of 
volts  per  meter,  while  the  function  0L(a;„2,  |/m2)  represents  the  phase  transformation 
of  the  lens  or  mirror  used  to  focus  the  light  collected  by  the  aperture  into  the  detector. 

Assuming  unity  pixel  hll-factor,  the  intensity  of  the  signal  at  a  pixel  of  the 
detector,  in  units  of  watts,  is  computed  by  forming  the  magnitude  squared  of  the 
held  propagated  to  the  detector, 


{^n3i  l/ms) 
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V^d^c\u^di^n3,ym3)\‘^, 


(2.5) 
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where  c  is  the  speed  of  light  in  a  vacuum,  e  is  the  permittivity  of  free  space,  and  is 
the  sampling  lattice  spacing  of  the  detector,  or  pixel  pitch.  The  detector  array  serves 
to  integrate  the  signal  over  some  discrete  number  of  coherence  times  r,  sample  the 
intensity  pattern  and  then  convert  the  signal  to  electrons. 

M 

yms)  ^  ,  (2-6) 

fc=l 

where  h  is  Plank’s  constant  and  M.  is  the  parameter  that  determines  the  degree  of 
temporal  and  spatial  speckle  averaging  that  occurs  due  to  the  limited  coherence  of 
the  laser  source  as  compared  to  the  duration  of  the  illumination  pulse. 

The  upper  limit  of  the  sum  in  Eqn.  2.6,  Af  is  a  parameter  that  indicates  the 
degree  of  coherence  of  the  optical  source  [31].  To  provide  range  gating  and  allow 
increased  signal-to-noise  ratio  of  the  received  illumination,  coherent  detectors  are 
often  gated  by  fairly  short  pulses.  Although  short,  the  length  of  the  gating  pulse  Tg 
is  often  many  time  longer  than  the  coherence  time  of  the  optical  source  such  that 
M.  =  Tgjr.  For  a  hxed  gate  time,  long  coherence  lasers  have  very  low  values  for 
Af,  while  narrowband  incoherent  light  sources  have  extremely  large  values  of  Af, 
with  commensurately  short  coherence  times,  r^.  A  more  complete  treatment  of  the 
speckle  parameter  includes  averaging  effects  due  to  spatial  correlation  in  addition  to 
purely  temporal  effects  as  described  above.  The  model  introduced  in  this  section 
uses  the  simplifying  assumption  that  spatial  correlation  effects  are  negligible.  The 
speckle  parameter  will  take  on  further  signihcance  in  Sec.  2.5  as  the  statistics  of  the 
illumination  are  considered  in  more  depth. 

2.2  Image  Sampling  for  Simulation  of  Coherent  Imaging 

The  creation  of  an  accurate  digital  representation  of  the  coherent  imaging  pro¬ 
cess  from  the  remote  scene  to  the  imaging  detector  requires  that  detailed  attention 
be  given  to  the  spatial  sampling  of  the  detector,  optical  aperture,  and  remote  target 
scene.  Additional  sampling  concerns  arise  when  simulating  accurate  statistical  phase 
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screens  necessary  to  model  the  tnrbulent  atmosphere,  as  will  be  further  discussed  in 
Sec.  2.3.  This  section  describes  the  framework  of  the  digital  simulation  that  is  nec¬ 
essary  to  model  the  experimentally  collected  imagery  described  in  Sec.  2.4.  Several 
example  calculations  are  included  that  use  the  actual  parameters  of  the  collection 
system  used  to  record  coherent  imagery  for  this  research  effort. 

Consider  the  simple  model  of  an  imaging  system  depicted  in  Fig.  2.1.  Reflected 
partially  coherent  light  from  the  target  with  a  mean  wavelength  of  A  is  propagated  over 
some  distance  L  to  the  optical  aperture  with  diameter  D.  The  light  is  subsequently 
propagated  to  the  imaging  detector  over  the  focal  distance  /  to  an  imaging  detector 
composed  of  a  discrete  array  of  detection  elements.  The  N  x  N  array  of  pixels  of  the 
imaging  detector  are  separated  by  a  spacing  of  in  both  axes. 

For  a  coherent  optical  system,  the  spatial  sampling  lattice  period  required  to 
satisfy  the  Nyquist  sampling  theorem  is  inversely  related  to  the  extent  of  the  optical 
aperture,  and  may  be  found  using  [30], 


Arf  = 


2D' 


(2.7) 


Using  example  dimensions  from  the  experimental  system  described  in  Sec.  2.4,  a  3 
meter  focal  length  and  20  cm  aperture  diameter  yield  a  minimum  detector  pixel  spac¬ 
ing  of  11.8  microns  given  laser  illumination  with  a  mean  wavelength  of  1.57  microns. 
It  should  be  noted  that  the  experimentally  collected  data  was  imaged  using  a  detector 
with  only  13  micron  pixel  spacing,  resulting  in  slightly  undersampled  imagery. 


Conversely,  the  Nyquist  required  sample  lattice  spacing  at  the  optical  aperture 
is  inversely  related  to  the  overall  extent  of  the  imaging  detector  based  on  Fresnel 
scaling  [30], 
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A/ 
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under  the  assumption  that  the  system  FOV  is  arranged  such  that  all  pixels  of  the 
detector  are  illuminated  by  the  remote  scene  image.  Using  the  hgures  above,  the 
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sample  spacing  at  the  aperture  is  approximately  0.7  millimeters,  requiring  at  least 
285  samples  per  axis  to  completely  sample  the  20  centimeter  aperture. 

A  wave-optics  simulation  was  created  to  adequately  model  the  propagation  of 
coherent  light  from  the  target  to  the  imaging  detector.  A  natural  choice  for  a  sampling 
lattice  at  the  aperture  is  to  match  the  physical  sampling  of  the  experimental  detector. 
Unfortunately,  using  the  entire  detector  grid  would  yield  unacceptably  high  computa¬ 
tional  requirements  for  generation  of  the  thousands  of  speckle  images  required  to  con¬ 
duct  turbulent  imagery  simulations,  given  the  complexity  of  the  Rayleigh-Sommerfeld 
propagation  integral  of  Eqn  2.3.  Reducing  the  held  of  view  by  a  factor  of  four  in  each 
direction  results  in  substantially  reduced  computational  complexity,  while  yielding 
an  optical  FOV  that  continues  to  dramatically  exceed  the  corresponding  isoplanatic 
patch  size  at  the  detector  for  all  but  the  very  weakest  of  turbulence  simulations.  This 
reduction  in  FOV  yielded  a  modest  128  x  128  detector  sampling  lattice  from  which 
simulation  parameters  were  derived  in  accordance  with  Fqn.  2.8. 

2.3  Atmospheric  Turbulence  Model 

Section  2.1  describes  a  model  that  is  well  suited  to  imaging  through  the  vac¬ 
uum  of  space  or  completely  undisturbed  air.  However,  any  helded  laser  vision  system 
would  require  optical  propagation  through  regions  of  atmosphere  corrupted  by  sig- 
nihcant  levels  of  turbulence.  To  simulate  the  effects  of  the  atmosphere,  one  well 
established  approach  is  the  treatment  of  the  continuous  atmospheric  path  between 
target  and  optical  aperture  as  a  series  of  discrete  thin  phase  screens  that  act  upon 
the  backscattered  coherent  illumination  from  the  target. 

A  key  research  goal  is  the  joint  estimation  of  the  scene,  together  with  param¬ 
eters  of  the  atmosphere  that  describe  the  resulting  image  blur.  The  estimation  of 
atmospheric  blur  parameters  may  be  considered  merely  as  a  by-product  of  the  scene 
estimation.  From  another  perspective,  atmospheric  condition  information  is  of  pri¬ 
mary  value,  as  it  allows  system  designers  and  operators  to  apply  this  information  to 
other  components  of  the  overall  system.  For  this  reason,  the  accurate  estimation  of 
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atmospheric  parameters  is  very  valuable.  Unfortunately,  image  data  collected  from  a 
remote  scene  is  not  often  accompanied  by  an  atmospheric  truth  source.  One  remedy 
is  to  collect  images  of  the  remote  scenes  in  concert  with  scintillometer  data  that  quan- 
tihes  the  atmospheric  turbulence.  In  practice,  such  measurements  are  difficult  and 
expensive  due  to  the  very  long  propagation  paths  typically  encountered.  Additionally, 
the  need  for  such  data  renders  previously  collected  scene  data  useless  for  atmospheric 
parameter  calibration  unless  other  means  are  employed  to  extract  such  information 
from  the  collected  image  data.  A  useful  atmospheric  parameter  estimation  technique 
that  avoids  scintillometry  measurements  is  offered  in  Sec.  2.9. 

The  accurate  construction  of  a  realistic  model  of  the  atmospheric  turbulence 
allows  parameterized  simulation  of  turbulence  degraded  imagery  that  can  be  processed 
by  the  joint  image/atmospheric  turbulence  parameter  estimator  algorithm  described 
in  Chap.  III.  Additional  insight  gained  from  the  atmospheric  model  may  be  used  to 
understand  the  relationship  of  the  phase-screen  correlation  and  the  effect  that  the 
atmosphere  has  on  causing  portions  of  the  observed  image  to  shift  spatially  at  the 
image  detector  plane.  The  correlation  between  these  shifts  is  an  important  tool  used 
in  the  derivation  of  an  anisoplanatic  optical  transfer  function,  as  outlined  in  Chap.  IV. 

It  is  assumed  that  the  variance  of  the  log-amplitude  fluctuations  at  the 
optical  aperture  is  sufficiently  small  that  the  effects  of  turbulence  are  dominated 
by  phase  effects  for  nominal  horizontal  path  imaging  scenarios.  This  assumption 
has  been  shown  to  be  useful  in  research  involving  atmospheric  turbulence  mitigation 
[61].  The  assumption  of  phase  dominated  atmospheric  conditions  is  required  for  the 
simplihcation  the  short  exposure  average  optical  transfer  functions  that  describes  the 
statistical  response  of  the  system  [31].  The  assumption  allows  a  farther  important 
simplihcation  that  eases  modeling  reqnirements.  By  ignoring  amplitude  scintillation 
effects  at  the  optical  aperture,  the  distinct  and  discrete  phase  screens  used  to  model  a 
volumetric  path  between  the  scene  and  lens  may  be  considered  as  a  single  thin  phase 
screen  placed  immediately  before  the  aperture.  Goodman  dehnes  the  random-phase 
screen  as  a  screen  that  “changes  the  phase  of  the  light  transmitted  in  an  unpredictable 
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Figure  2.2:  Phase  screen  representation  of  turbulent  atmosphere.  The  volume  of 

turbulent  atmosphere  is  broken  into  a  hnite  number  of  discrete  phase  volumes.  Each 
of  these  component  volumes  is  represented  by  a  distinct  phase  screen.  The  screens 
may  be  collapsed  to  a  single  phase  screen  located  at  the  aperture  given  assumptions 
regarding  the  strength  of  the  turbulence. 

fashion  but  does  not  appreciably  absorb  light”  [31].  The  summed  contribution  of  the 
individual  phase  screens  must  be  carefully  considered  to  allow  creation  of  a  physically 
accurate  model. 

Figure  2.2  depicts  the  treatment  of  continuous  phase  perturbations  caused  by  the 
atmosphere  as  a  hnite  number  of  discrete  phase  screens,  each  with  spatial  correlation 
properties  that  are  described  by  atmospheric  conditions.  Such  a  simplihed  model 
allows  the  generation  of  arbitrary  images  through  various  levels  of  random  turbulence. 

To  construct  the  model,  the  characteristics  of  the  turbulence  must  be  repre¬ 
sented  by  the  individual  phase  screens.  Turbulence  is  most  conveniently  modeled  by 
quantifying  the  statistical  distribution  of  the  random  turbulent  eddies  caused  by  the 
distributed  heating  and  cooling  of  the  atmosphere.  The  power  spectral  density  (PSD) 
of  the  phase  screen  may  be  quantihed  using  several  useful  atmospheric  tur¬ 

bulence  models,  including  the  Von  Karman  PSD  [60] 


(2.9) 
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where  Kq  =  2ti/ Lq  and  Km  =  5.92//o  with  Lq  representing  the  enter  scale  diameter  of 
the  tnrbnlence,  and  /q  the  inner  scale  diameter,  in  meters.  C‘^{z)  is  the  atmospheric 
strnctnre  constant  and  represents  the  mnltiplier  of  the  PSD  to  acconnt  for  the  overall 
strength  of  the  tnrbnlence.  Typical  valnes  for  C‘^{z)  might  range  from  10“^^ 
nnder  excellent  nighttime  seeing  conditions  to  nnder  extremely  poor 

daytime  atmospheric  tnrbnlence  conditions  [61].  The  fnnctional  dependence  of  C^^z) 
on  distance  from  the  apertnre  may  be  dropped  if  a  homogeneons  tnrbnlence  volnme 
is  assnmed,  as  might  be  enconntered  for  horizontal  imaging.  In  this  case, 
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Essentially,  the  physical  model  approximates  the  optical  delay  contribnted  by  a 
single  layer  of  the  atmosphere  as  a  Ganssian  random  process  with  an  antocorrelation 
fnnetion  that  depends  only  on  differences  between  locations  in  the  apertnre  of  the 
telescope  [31].  Following  Roggemann  and  Welsh  [60],  if  (xi,  yi)  and  {x2, 1/2)  are  spatial 
locations  in  a  plane  containing  a  random  phase-screen,  then  the  relative  correlation 
distance  between  two  parts  of  the  phase-screen  is  dehned  by  pm  =  [\x2  —  xi\‘^  + 
\y2  —  1/1 Let  (f)m  be  the  optical  delay  in  the  pnpil  plane  in  radians  dne  to  the 
layer  of  the  atmosphere  in  nnits  of  waves.  A  key  qnantity  in  characterizing  the 
statistics  of  the  phase  screen  is  the  antocorrelation  of  0^,  denoted  by  Rm-  Under  the 
assnmption  of  atmospheric  isotropy  and  homogeneity,  the  fnnetion  Rm  depends  only 
on  the  relative  radial  distance,  pm,  between  two  locations  of  the  phase-screen,  i.e.. 


Rm{Xl,yi,X2,y2)  =  E[(l)m{xi,yi)(j)m{x2,y2)]  =  Rm{Pm)-  (2.11) 


The  single-layer  model  described  above  can  be  extended  to  describe  the  phase 
distortion  indneed  by  the  entire  atmosphere  throngh  a  mnlti-layer  model  [60].  Be- 
canse  the  atmospheric  layers  are  assnmed  to  be  statistically  independent,  the  cross¬ 
correlation  of  the  optical  delay  at  two  points  in  time  and  space  is  the  snm  of  the 
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cross-correlation  of  each  of  the  layers  [60], 

M 

yii  X2i  y2)  =  ^  ^  RmjPm}-:  (2-12) 

m=l 

where  M  is  the  number  of  layers.  If  the  von  Karman  power  spectral  density  of  Eq.  2.10 
is  used  to  describe  the  spatial  distribution  of  the  index  of  refraction  in  each  of  the 
2-D  phase  screens,  the  expression  for  the  autocorrelation  using  an  M-layer  model  is 
then  given  by  [60], 

,2.13) 

m=l 

where  is  a  modihed  Bessel  function  of  the  second  kind  of  order  5/6,  A  is  the  wave¬ 
length  of  the  light,  5zm  is  the  thickness  of  each  atmospheric  volume  in  the  direction 
of  propagation,  and  Lo  is  the  outer  scale  of  the  turbulence. 

The  size  of  the  random  phase  screen  is  determined  by  the  size  of  the  column  of 
atmosphere  originating  at  each  point  on  the  target  and  propagating  to  the  telescope 
aperture.  This  is  determined  by  the  size  of  the  aperture,  the  location  of  the  extreme 
points  of  the  scene  under  observation  and  the  spatial  sample  rate  in  the  aperture  of 
the  telescope.  Additionally,  care  must  be  taken  to  properly  size  the  overall  extent 
of  the  phase  screens  due  to  sampling  considerations.  The  lowest  spatial  frequency 
component  attainable  in  simulation  is  inversely  proportional  to  the  overall  extent  of 
the  phase  screen.  In  order  to  properly  account  for  low  frequency  contributions  of 
atmospheric  turbulence  such  as  tip  and  tilt,  relatively  large  phase  screens  must  be 
constructed  at  the  expense  of  large  computational  burden  [35]. 

This  autocorrelation  model  can  be  used  to  generate  random-phase  screens  that 
possess  the  proper  spatial  correlations.  Given  the  appropriate  spatial  autocorrelation 
function  for  each  layer,  individual  realizations  of  turbulence  may  be  generated  by 
producing  a  matrix  of  zero  mean,  unit  variance  Gaussian  random  numbers,  G{xi,  yi). 
Random  2-D  uncorrelated  phase  screens  may  be  appropriately  hltered  using  the  ap- 
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propriate  correlation  function  by  pointwise  multiplying  the  Fourier  transform  of  the 
the  derived  autocorrelation  function  R{xi,yi,X2,y2)  with  the  spectral  domain  Gaus¬ 
sian  screen  [38]. 

A  single  realization  of  the  turbulence  for  a  single  layer  may  be  produced  from 
the  inverse  2-D  Fourier  transform  of  the  multiplication  of  Fg{u,v),  which  is  the  2-D 
Fourier  transform  of  the  zero- mean,  unit  variance  Gaussian  random  matrix  G{xi,  yi), 
and  Fm{u,v)  which  is  the  square  root  of  the  2-D  Fourier  transform  of  the  correlation 
function  Rm{Pm), 

Ni  Ni 

M^,y)  =  J^Y.FG{u,v)FUu,v)e-^^-^^^+^y'>/^K  (2.14) 

u=l  V  =  1 

This  process  yields  a  phase  screen  of  Nf  number  of  points  that  statistically  possesses  a 
spatial  autocorrelation  equal  to  R^,  as  long  as  the  matrix  containing  R^  is  constructed 
to  be  larger  than  twice  the  outer-scale  of  the  turbulence,  Lq,  and  is  sampled  at 
less  than  the  period  of  the  inner-scale  of  the  turbulence,  /q.  The  total  phase  screen 
(l>a{xn,  ym,  Xn2,  l/m2)  between  a  point  in  the  target  plane  (a;„,  yn)  and  a  point  {xn2,  ym2) 
in  the  aperture  plane  may  be  computed  by  summing  the  contributions  along  the 
unique  path  through  individual  phase  screens  in  each  layer  between  the  target  and  the 
aperture.  This  path-dependent  summation  necessitates  the  creation  of  a  new  phase 
screen  for  each  point  {xn,yn)  propagated  from  the  target  to  the  optical  aperture.  The 
hnal  result  is  a  collection  of  random-phase  screens  that  encodes  the  proper  degree  of 
spatial  correlation  based  on  the  strength  and  scale  of  the  atmospheric  turbulence,  for 
a  given  optical  path  from  each  target  point  source  to  the  optical  aperture. 

The  creation  of  an  aggregate  phase  screen  allows  the  inclusion  of  stochastic 
turbulence  effects  by  multiplication  of  the  field  arriving  at  the  optical  aperture  given 
in  Eqn.  2.3  with  a  thin  lens  transmissivity  function  described  by  the  phase  relationship 
of  Eqn.  ??.  The  field  propagated  to  the  imaging  detector  plane  becomes 
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where  Tiens  =  and  Tgcreen  = 

Reference  to  Fig.  2.2  demonstrates  that  the  aggregate  summed  random-phase 
screen  developed  to  model  the  imaging  of  point  A  will  be  considerably  different  than 
that  employed  to  image  point  B  or  C.  This  is  expected  and  desired  for  a  model  of  a 
system  that  exceeds  the  isoplanatic  angle  described  in  Sec.  1.3.  Rayleigh-Sommerfeld 
propagation,  while  computationally  expensive,  preserves  the  path-dependent  rela¬ 
tionships  necessary  to  simulate  anisoplanatic  turbulence  effects  for  a  spatially  variant 
imaging  system. 


2.4  Experimental  Data  Collection  System 

An  experimental  coherent  imaging  system  was  assembled  by  the  Air  Force  Re¬ 
search  Laboratories  (AFRL)  Sensors  Division.  Recent  advances  in  the  design  of 
electron-bombarded  charge-coupled  device  (EBCCD)  imaging  microchips  has  enabled 
the  efficient  capture  of  photons  in  the  near  infrared  region  used  by  the  imaging  system. 
The  brassboard  system  is  shown  in  Fig.  2.3. 

The  coherent  imaging  system  was  used  to  collect  a  very  extensive  set  of  target 
images  at  a  variety  of  ranges  between  3  and  27  kilometers.  Table  2.1  describes  some 
important  operating  characteristics  of  the  laser  used  to  illuminate  the  remote  scene, 
while  Table  2.2  describes  the  optical  receiver  system  [47]. 

The  individual  images  obtained  by  the  coherent  optical  system  were  heavily 
corrupted  by  the  speckle  that  is  caused  by  the  random  variations  of  surface  roughness 
that  are  on  the  order  of  an  optical  wavelength  or  larger,  as  well  as  atmospheric 
turbulence  speckle.  Figure  2.4  shows  a  representative  speckle  image  of  a  target  imaged 
at  a  range  of  10  kilometers.  To  combat  the  effects  of  image  speckle,  a  series  of 
successive  images  were  registered  (motion  compensated)  and  averaged.  Figure  2.5 
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Figure  2.3:  Brassboard  coherent  imaging  system.  A  pulsed  laser  illuminator  with 
relatively  long  coherence  time  is  aimed  along  the  optical  boresight  of  a  small  aperture 
Cassegrain  telescope.  The  system  was  used  to  collect  a  large  database  of  coherently 
illuminated  images  of  various  targets  surrounding  a  mountaintop  test  range. 

shows  a  representative  image  formed  by  registering  a  series  of  coherent  speckle  images 
collected  by  the  laser  vision  system.  For  imagery  collected  at  a  range  of  approximately 
10  kilometers,  the  system  has  an  instantaneous  £eld-of-view  (FOV)  of  approximately 
2.23  milli-radians  (using  a  three  meter  focal  length  lens).  This  wide  FOV  is  several 
orders  of  magnitude  larger  than  the  isoplanatic  angle  for  a  typical  atmospheric  profile, 
as  calculated  using  expressions  presented  in  [60].  This  suggests  that  the  isoplanatic 
patch  size  of  the  image  detector  array  ranges  from  several  tens  of  pixels  to  as  small 
as  a  single  pixel  or  less.  It  is  also  assumed  that  the  turbulence  seen  by  each  frame 
acquired  by  the  laser  vision  system  will  be  statistically  independent.  This  assumption 
is  based  on  Taylor’s  Frozen  Flow  hypothesis  [60]  in  conjunction  with  the  relatively 
small  aperture  diameter  of  the  system,  the  10  Hz  frame  rate  and  surface  winds  usually 
in  excess  of  10  Knots  during  the  testing  of  the  sensor. 

The  required  inter-frame  period  discussed  above  strictly  applies  to  each  inter¬ 
mediate  thin  phase-screen  representations  along  the  path  to  the  target.  However  the 
atmosphere  closest  to  the  aperture  creates  the  most  significant  phase  aberrations  at 


2-15 


Table  2.1;  This  table  describes  the  parameters  of  the  laser  system  used  to  illuminate 
remote  targets  at  the  North  Oscura  Peak  (NOP)  site  of  the  White  Sands  Missile 
Range,  New  Mexico. 


Parameter 

Value 

Laser  Type 

Nd:YAG  w/KTP  Optical  Parametric  Oscillator 

Wavelength 

1.57  micrometers 

Pulsewidth 

12  nanoseconds 

Pulse  Repetition  Rate 

up  to  20  Hz  (10  Hz  limit  due  to  camera) 

Output  Energy 

130  milli Joule/pulse  at  1.57  micrometers 

Beam  Diameter  at  1/e 

13  millimeters 

Beam  Divergence  at  1/e 

15  milliradian  at  1.57  micrometers 

the  imaging  plane.  For  the  case  of  an  airborne  sensor,  this  limitation  is  not  severe, 
as  simple  calculations  reveal  that  the  aperture  movement  through  the  atmosphere 
need  not  be  great  for  nominal  sizes  of  the  outer  scale  of  turbulence,  Lq.  However, 
on  a  ground-based  optical  range,  wind  speeds  should  be  sufficiently  high  to  satisfy 
this  requirement.  For  example,  for  the  case  of  a  1-meter  outer  scale  turbulence  size, 
a  transverse  wind  speed  of  10  m/s  would  result  in  completely  refreshed  (and  thus 
statistically  uncorrelated)  atmosphere  preceding  the  aperture  every  10  ms.  A  frame 
rate  of  10  Hz  would  admit  this  condition.  In  the  case  of  the  experimental  collection 
system,  the  gated  laser  detector  charge-coupled  device  (CCD)  has  a  very  fast  shutter 
speed  of  approximately  12  ns  and  a  frame  capture  rate  of  10  Hz.  Over  the  course 
of  the  optical  testing  for  the  data  presented  in  Chapters  HI,  IV  and  V,  wind  speed 
remained  above  35  m/s  with  gusts  much  higher  than  that  hgure  over  the  duration 
of  the  test  period.  A  10  Hz  frame  rate  would  thus  permit  outer  scales  of  turbulence 
on  the  order  of  3.5  meters.  Such  a  hgure  provides  a  nominal  value  of  expected  tur¬ 
bulence  outer  scale  according  to  the  literature  [60].  Although  there  may  be  a  small 
degree  of  correlation  between  screens  generated  for  different  frames,  the  assumption 
of  independence  greatly  simplihes  the  turbulence  simulation. 
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Figure  2.4:  Sample  image  collected  by  the  coherent  vision  system.  The  image 

suffers  from  heavy  degradation  due  to  coherent  speckle  effects  caused  by  target  surface 
roughness. 
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Figure  2.5:  Averaged  image  collected  by  the  coherent  vision  system.  The  image 

was  created  by  averaging  50  spatially  registered  frames.  Registration  and  averaging 
is  typically  required  to  remove  the  objectionable  effects  of  coherently  formed  speckle 
that  is  present  in  each  of  the  individual  frames. 
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Table  2.2;  This  table  describes  the  parameters  of  the  laser  vision  receiver  imaging 
system  used  to  collect  experimental  imagery  of  remote  targets  at  the  North  Oscura 
Peak  site  of  the  White  Sands  Missile  Range,  New  Mexico. 


Parameter 

Value 

Imaging  Detector 

Electron  Bombarded  Charge  Coupled  Device 

CCD  Manufacturer 

Intevac,  Inc. 

CCD  Wavelength 

1.57  micrometers 

CCD  Quantum  Efficiency 

approximately  30% 

CCD  Output  Data 

12-bit  digital 

CCD  Array  Size 

512x512  pixels 

CCD  Pixel  Spacing 

13  X  13  micrometers 

CCD  Active  Imaging  Area 

6.7  X  6.7  millimeters 

Optical  Filter 

1.57  micrometer  laser  line  hlter 

Optical  Aperture 

8  inch 

Optical  Telescope  Manufacturer 

Celestron 

Focal  Length 

2000  millimeter  with  1.5x  and  2x  extenders 

F-number 

9.84  basic,  14.76,  19.68,  29.52  with  extenders 

Frame  Capture  Rate 

10  frames  per  second 

2.5  Image  Speckle  Parameter  Estimation 


One  of  the  requirements  of  scene  and  atmospheric  parameter  estimation  is  to 
hrst  obtain  an  estimate  of  the  degree  of  temporal  coherence  of  the  laser  source  used 
to  illuminate  the  remote  scene.  The  following  section  documents  a  novel  estimator 
developed  over  the  course  of  this  research  effort  to  estimate  the  speckle  parameter,  Ai 
from  a  series  of  collected  images  taken  at  an  experimental  optical  range.  The  speckle 
parameter  is  essentially  a  free  parameter  that  must  be  estimated  prior  to  employing 
the  image  restoration  techniques  described  in  Chapters  I  and  IV. 


The  photon  distribution  for  individual  pixels  of  images  collected  by  partially 
coherent  imaging  systems  has  been  shown  to  follow  a  negative  binomial  random  pro¬ 
cess  [31]. 


P{K) 


T{K  +  M) 
T{K  +  1)T{M) 


1  +  ^' 
K 


-K  r 


K 


-M 


(2.16) 


2-19 


0.05 


Figure  2.6:  Plot  of  the  negative  binomial  probability  mass  function  with  Kavg  =  50. 
Curves  are  plotted  for  values  of  M.  ranging  from  2  to  100.  Note  that  in  the  limit 
as  A4  increases  without  bound,  the  PMF  approaches  that  of  the  familiar  Poisson 
distribution. 


where  the  distribution  is  parameterized  by  the  mean  photon  count  collected  within  a 
particular  pixel,  K,  and  the  speckle  parameter  Ai.  The  speckle  parameter  may  also 
be  regarded  as  an  indication  of  the  variance  of  the  distribution.  Figure  2.6  shows 
various  probability  mass  functions  for  differing  values  of  the  speckle  parameter,  given 
a  constant  mean  of  =  50  photons.  In  the  limit  as  At  grows  without  bound, 
the  variance  decreases  to  approach  the  mean  of  the  distribution.  This  behavior  is 
expected,  as  the  negative  binomial  distribution  tends  to  Poisson  as  Ai  grows  large. 

The  derivation  of  a  maximum  likelihood  estimator  for  the  speckle  parameter  Ai 
is  straightforward  [71].  The  following  development  assumes  integer  values  of  Ai  to 
expedite  calculation  of  a  useful  estimate.  For  the  case  of  J  observations  of  a  single 
pixel  of  an  image,  one  can  write  the  distribution  as  conditional  on  the  two  non-random 
parameters,  K  and  Ai 


P{K  I  K,Ai)=ll 


T{Kj +  Ai)  r  At 
r  {Kj  +  1)  F  (M)  f  ^  K 


-Ki 


K 


H  -M 


(2.17) 
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For  convenience,  the  natural  logarithm  of  the  conditional  probability  may  be 


computed  as 


In  [P(K  I  K,M)]=Y. 

J  =  1 


V{K^  +  M) 

r  (iF,- + 1)  r  (Af ) 


iF,ln 


-Alin 


(2.18) 


The  derivative  of  Eqn.  2.18  may  be  easily  shown  to  be 


d 

dK 


{ln[P(A' |/?,VI)]}=5;|a' 

i=i 


M 

MK  +  K^ 


M  \ 
K  +  M  j  ’ 


(2.19) 


Setting  Eqn.  2.19  to  zero  and  solving  for  a  zero  yields  a  candidate  for  the  maximum 
likelihood  estimate  for  the  pixel  intensity 


1  ^ 

^rnl  =  -j  Kj . 

i=i 


(2.20) 


It  is  easy  to  show  that  K^i  yields  a  global  maximum  for  the  likelihood.  The  estimate 
is  not  surprising  and  may  easily  be  shown  to  be  unbiased.  A  more  useful  estimator 
may  be  found  for  the  speckle  parameter  M..  However,  the  maximization  approach 
used  above  is  intractable,  and  a  numerical  maximization  of  the  log-likelihood  equation 
of  Eqn  2.18  is  practical. 

Since  T  (n  -|-  1)  =  n!  for  positive  integer  n,  it  follows  that 

k-l 

lnr(A;)  =  ln(A;- 1)!  =  ^Ini  (2.21) 

i=\ 
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valid  for  all  k  >  1.  Using  this  relationship,  Eqn.  2.18  may  be  rewritten  as 


In  [P  {K  I  K,M)] 


J  (  Kj+M-l 


Ki 


M-l 


ln*-^ln*-  ^In* 

j=i  I  i=i 


2=1  2  =  1 


-  K,  In 


1  +  ^' 
K 


-Min 


K' 


(2.22) 


The  computational  burden  of  Eqn.  2.22  is  reduced  by  the  observation  that  the  hrst 
Kj  summation  terms  of  the  hrst  two  terms  within  brackets  cancel,  leaving 


In  [P  {K  I  K,M)] 


j  (  Kj+M-l  M-l 

j=l  f  i=Kj+l  i=l 


Ini  —  Ini 


-  K,  In 


-Min 


(2.23) 


The  closed  form  estimate  of  the  mean  intensity  found  in  Eqn.  2.20  may  be  substituted 
into  Eqn.  2.23  and  the  result  maximized  using  standard  numerical  techniques  such  as 
steepest  decent  or  Newton’s  method. 

The  utility  of  the  estimator  described  above  lies  in  its  application  to  experimen¬ 
tally  collected  imagery.  For  a  laser  vision  system  composed  of  an  illuminator  with 
hxed  integration  period  or  gating  time  Tg  and  stable  coherence  time  r,  the  speckle 
parameter  may  be  considered  unknown  but  constant  over  reasonable  collection  pe¬ 
riods,  given  stable  operating  temperatures  and  other  system  factors.  Given  a  large 
ensemble  of  experimentally  collected  image  data,  a  useful  estimate  of  the  system  M 
may  be  calculated  using  the  estimator  above. 

Laser  illumination  propagated  from  the  platform  emitter  to  a  remote  scene  will 
undergo  some  degree  of  amplitude  and  phase  variation  due  to  atmospheric  turbulence 
effects  discussed  in  Sec.  2.1.  It  is  assumed  that  this  additional  variance  produces  an 
image  ensemble  that  is  also  modeled  by  the  negative  binomial  distribution,  despite 
amplitude  correlation  effects  observed  at  the  target  due  to  the  size  of  the  laser  aperture 
or  the  turbulent  seeing  conditions. 
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2.6  Effect  of  Atmospheric  Amplitude  and  Phase  Distortion  on  Intensity 
Variation  at  the  Detector 

In  the  vacuum  of  space,  it  can  be  mathematically  demonstrated  that  the  in¬ 
tensity  variation  of  the  detected  image  formed  under  partially  coherent  illumination 
closely  follows  a  negative  binomial  stochastic  distribution  [31].  A  simple  model  of  the 
coherent  illumination  of  a  remote  target  by  a  laser  beam  was  presented  in  Sec.  2.1.  It 
is  not  clear,  however,  that  this  relationship  holds  for  the  case  of  coherent  propagation 
through  a  large  volume  of  random  turbulence,  as  might  be  encountered  in  terrestrial 
optical  imaging  applications.  The  following  analysis  will  demonstrate  why  such  a  sim¬ 
ple  model  is  adequate  for  the  purposes  of  this  research  study.  The  analysis  will  begin 
with  a  discussion  of  random  phasor  sums  as  required  to  understand  the  amplitude 
and  phase  fluctuations  of  a  held  propagated  through  random  atmospheric  turbulence. 
The  discussion  will  then  be  extended  to  the  negative  binomial  distribution  discussed 
in  Sec.  2.5  and  it  will  be  shown  that  the  resulting  distribution  arriving  at  the  detector 
remains  negative  binomial,  despite  random  amplitude  and  phase  huctuations  that 
affect  the  coherent  illumination  rehected  from  the  target. 


2. 6. 1  Random  Phasor  Sums.  The  following  analysis  describes  the  statistics 
of  the  amplitude  and  phase  of  a  large  sum  of  random  phasors,  and  closely  follows 
the  development  of  [18].  The  discussion  is  motivated  by  the  Rayleigh-Sommerfeld 
integrals  used  for  propagation  of  coherent  light  from  the  target  to  the  aperture  and 
detector.  For  example,  referring  to  Equations  2.2  and  2.3,  it  is  clear  that  each  point  on 
the  optical  aperture  is  formed  by  the  summation  of  many  such  phasors,  each  of  which 
has  a  random  amplitude  and  phase  due  to  target  roughness  under  typical  optical 
imaging  conditions.  Given  a  large  set  of  N  vectors  having  random  amplitude  and 
random  phase  (j)k,  one  may  construct  the  normalized  random  phasor  sum 


N 


a  =  = 


.j4>k 


(2.24) 
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from  which  the  real  and  imaginary  components  may  be  expressed  as 


r  =  3?{a} 
i  =  9{a} 


1  .y  ^ 

1  ^ 

—  ^\ak\sin{<i)k). 

'  k=l 


(2.25) 


The  expected  value  of  each  component  is  easily  calculated  since  the  expectation 
of  both  the  sin  and  cos  function  is  zero  for  uniform  phase, 


(^)  ^(|afc|)(cos(0fc))  =  0 

1  ^ 

(i)  ^(|afc|)(sm(0fc))  =  0.  (2.26) 

Furthermore,  the  real  and  imaginary  parts  of  the  random  phasors  are  uncorre¬ 
lated,  since  {cos{(f)k) sin{(f)rn))  =  0, 

^  N  N 

{\ak\\am\){cos{(j)k)sin{(j)m))  =  0.  (2.27) 

k=l  m=l 

The  variance  of  the  real  and  imaginary  components  is  identical. 


^  N  N 

(|afc||am|)(cos(0fc)cos(0m)) 

k=l  m=l 
AT  TV 

( Ittfc  1 1  a™  I )  {sin{(t)k)sin{(t)m)) 

k=l  m=l 


1  ^  (b.P) 

N  ^  2 

k=l 


N  ^  2 

k=l 


(2.28) 


where  {\oik\^)  is  the  second  moment  of  the  amplitude  of  the  phasors  and  is  a  property 
of  the  target  reflectance. 
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Given  a  large  number  of  random  phasors  where  the  real  and  imaginary  compo¬ 
nents  are  independent,  the  Central  Limit  Theorem  allows  construction  of  a  joint  PDF 
of  the  resulting  summation  which  is  circularly  Gaussian 

where 

1  ^  /I  I 

=  lim  (r^)  =  lim  =  lim  —  — .  (2.30) 

Af^oo'  '  V^oo'  '  AT^oo  N  ^  2  ^ 

k=l 

This  result  is  strictly  applicable  only  in  situations  where  the  random  amplitudes 
ak  are  independent  and  the  phase  <pk  is  uniformly  distributed  on  the  interval  [— 7r,7r]. 
Although  it  might  seem  an  improbable  distribution,  the  assumption  of  uniform  phase 
is  appropriate  to  the  imaging  physics.  For  an  unknown  surface  roughness  of  a  target 
that  has  a  variance  several  times  larger  than  the  mean  optical  wavelength  of  the  illu¬ 
minating  coherent  beam,  the  resulting  random  phase  becomes  uniformly  distributed 
on  the  interval  [— vr,  tt]  due  to  phase  wrapping  effects. 

Since  the  amplitude  and  phase  are  related  to  the  real  and  imaginary  components 
of  the  phasor,  r  and  i  by  the  relationships 


a  = 

9  =  tan“^  -,  (2.31) 

r 


the  joint  density  function  of  a  and  6  may  be  found  using  the  Jacobian  matrix  and 
may  be  expressed  as 


PAe{a,0) 


I  ^expj-^} 
0 


— TT  <  6  <  TT, 

a  >  0, 

otherwise. 


(2.32) 
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from  which  the  marginal  distributions  of  a  and  6  may  be  easily  calculated  by  integra¬ 
tion  [31] 


PA{a) 


exp 


0 


a  >  0, 
otherwise, 


Pe{0) 


^  -n<e<7i, 

0  otherwise. 


(2.33) 


The  marginal  densities  of  a  and  6  in  Eqn.  2.33  may  be  multiplied  to  obtain  the 
joint  density  of  Eqn.  2.32,  which  demonstrates  that  A  and  0  are  also  independent 
random  variables  as  are  the  real  and  imaginary  components  of  the  phasors,  as  noted 
above  by  Eqn.  2.27. 

The  analysis  above  is  central  to  the  argument  that  the  distribution  of  intensity 
at  the  detector  plane  follows  a  negative  binomial  distribution  despite  the  effects  of 
turbulence.  To  understand  this  statement,  it  is  first  necessary  to  examine  several  of 
the  assumptions  that  lead  to  a  negative  binomial  model  for  the  detected  intensity.  For 
an  integration  time  Tg  much  shorter  than  the  coherence  time  r  of  the  optical  illumina¬ 
tion,  a  detected  pixel  intensity  is  formed  by  a  random  phasor  sum  as  discussed  above, 
and  has  a  complex  envelope  amplitude  A  which  is  Rayleigh  as  shown  in  Eqn.  2.33. 
Since  the  intensity  I  is  the  square  of  the  amplitude,  the  distribution  transforms  to 
negative  exponential  [31] 


Pi{I)  = 


(2.34) 


for  /  >  0  and  zero  otherwise.  However,  for  an  arbitrary  counting  interval  much 
longer  than  the  laser  coherence  time,  many  such  negative  exponential  distributions 
occur  during  the  integration  period  of  the  detector,  Tg.  Defining  W  as  the  integrated 
intensity  over  the  counting  period,  the  distribution  of  W  can  be  shown  to  follow  a 
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gamma  probability  density  function  [31], 


Pw{W) 


(2.35) 


where  W  is  the  mean  integrated  intensity,  and  M.  is  the  degree-of-freedom  of  the  dis¬ 
tribution,  identical  to  the  speckle  parameter  estimated  in  Sec.  2.5.  The  detection  of 
intensity  follows  a  doubly-stochastic  process  characterized  by  the  statistics  of  the  illu¬ 
mination  Pw(W),  as  well  as  the  statistics  of  the  photon-matter  interaction,  P{K\W). 
The  latter  is  commonly  modeled  to  follow  a  Poisson  distribution.  Mandel’s  formula 
may  be  used  to  find  the  unconditional  photon  distribution  at  the  detector  plane  [31], 
where 

OO 

p{K)  =  I  P{K\W)pw{W)dW 
b 

OO 

=  j  ^^^e-^^pw{W)dW,  (2.36) 

0 


and  a  is  related  to  the  mean  integrated  intensity  hy  K  =  aW. 

Substitution  of  Eqn.  2.35  into  Mandel’s  formula  of  Eqn.  2.36  yields  the  negative 
binomial  distribution  of  intensity  at  the  detector  for  an  arbitrary  counting  interval 
given  in  Eqn.  2.16  and  repeated  below  for  convenience. 


P(K) 


r(K  +  M) 
r(K  +  i)r(M) 


1  +  ^' 
K 


-K 


K 


H  -M 


The  key  assumptions  used  to  develop  the  random  phasor  sum  analysis  were 
independence  of  the  amplitudes  and  phases  of  the  individual  phasors  to  be  summed. 
The  fact  that  the  real  and  imaginary  parts  of  the  resulting  amplitude  phasors  at  the 
optical  aperture  are  uncorrelated,  together  with  the  result  that  the  phase  is  uniformly 
distributed  on  the  interval  [— 7r,7r),  leads  to  the  important  observation  that  the  same 
conditions  are  satished  for  the  subsequent  propagation  from  aperture  to  detector.  The 
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introduction  of  a  random  phase  due  to  atmospheric  turbulence  does  not  affect  this 
analysis,  since  the  phase  due  to  the  random  phasor  sum  is  uniform  at  the  aperture. 
The  addition  of  further  phase  effects  will  wrap  modulo  27r  and  result  in  uniform  phase 
distribution  at  the  aperture.  The  result  is  a  negative  binomial  intensity  distribution 
at  the  detector,  despite  the  addition  of  varying  degrees  of  random  phase  disturbance 
due  to  atmospheric  turbulence. 


2. 7  Effect  of  Image  Intensity  Scaling  and  Quantization 

The  estimator  development  of  Chap.  Ill  is  based  on  the  statistics  of  photon 
arrival  at  the  imaging  detector  plane.  All  digital  image  collection  systems  impose 
distortion  on  the  statistics  of  the  captured  images  due  to  the  unavoidable  scaling 
and  quantization  of  the  analog-to-digital  (A/D)  process,  unless  the  added  expense  of 
a  true  photon-counting  detector  is  justihed.  This  distortion  may  be  regarded  as  a 
modification  of  the  probability  distribution  of  the  intensity. 

Pixel  intensity  scaling  that  often  accompanies  the  quantization  process  must  be 
removed  by  characterization  of  the  optical  system.  In  most  systems,  this  entails  scal¬ 
ing  of  pixel  intensity  by  a  factor  that  represents  the  number  of  photons  per  intensity 
count.  This  factor  is  unity  in  a  photon-counting  camera  system,  but  may  be  quite 
high  for  less  sensitive  detectors.  Without  scale  correction,  the  intensity  statistics  of 
partially  coherent  illumination  can  no  longer  be  accurately  modeled  by  a  negative 
binomial  random  process.  Consider  a  scaling  factor  of  g  =  1/p,  where  p  represents 
the  number  of  photons  per  intensity  count  stored  by  the  detector  system.  The  mean 
of  the  photon  intensity  distribution  is  K,  while  the  variance  is  [31] 


K 
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(2.37) 


After  scaling,  the  detected  mean  becomes  qK,  while  the  detected  variance  is 
which  is  not  equal  to 

=  qK{l+'i^y  (2.38) 

as  required  by  specihcation  of  the  negative  binomial  process.  For  non-photon  counting 
detectors,  g  <  1,  and  the  detected  variance  is  lower  than  that  predicted  of  a  true 
negative  binomial  process.  Unless  properly  compensated  for,  this  phenomenon  yields 
erroneously  high  estimates  of  the  speckle  parameter  from  partially  coherent  image 
frames. 

As  will  be  discussed  in  Sec.  3.2.3,  accurate  formulation  of  the  statistical  process 
variance  is  necessary  to  conclude  iterations  of  the  scene  estimation  algorithm.  For 
non  photon-counting  systems  the  scaled  image  variance  is  calculated  by  q^a]^  using 
Eqn.  2.37. 

Even  after  device  scaling  calibration  is  completed,  the  effect  of  pixel  intensity 
quantization  remains.  In  systems  with  coarse  A/D  intensity  quantization,  relatively 
large  intensity  variance  is  masked  by  the  coarse  assignment  of  these  fluctuations  to 
relatively  few  intensity  bins.  This  phenomenon  is  especially  notable  in  darker  pixels 
where  image  SNR  is  low  and  the  standard  deviation  of  the  pixel  intensity  is  small 
compared  to  the  quantization  effect.  The  apparent  variance  is  limited  with  a  corre¬ 
sponding  artihcial  increase  in  the  estimated  speckle  parameter  per  Eqn.  2.37.  There¬ 
fore,  estimation  of  the  laser  illuminator  speckle  parameter  as  discussed  in  Sec.  2.5  is 
best  accomplished  using  bright  pixels  with  relatively  large  means.  A  bright  pixel  is 
therefore  dehned  as  one  whose  mean  is  high  enough  that  the  standard  deviation  of 
intensity  is  at  least  twice  the  quantization  step  size.  Such  a  condition  ensures  that 
the  noise  process  is  dominated  by  intensity  variation  versus  quantization  effects.  Fur¬ 
thermore,  darker  pixels  from  experimentally  collected  images  tend  to  be  governed  by 
random  processes  such  as  additive  noise  introduced  by  system  amplihers  or  shot-noise 
typical  of  solid-state  detector  systems. 
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2.8  Effects  of  Image  Registration  and  Averaging 

The  relatively  long  coherence  times  of  candidate  laser  illuminators  cause  a  signif¬ 
icant  degree  of  laser  speckle  for  individually  collected  frames  of  2-D  LIDAR  imagery. 
Additionally,  the  instantaneous  distortion  of  the  turbulent  atmosphere  causes  what 
can  be  described  as  atmospheric  speckle.  Figure  2.7(a)  illustrates  these  phenomena 
given  a  representative  imaging  scenario.  Although  created  by  distinct  random  pro¬ 
cesses,  the  cumulative  effect  of  these  distortions  may  be  effectively  reduced  by  simple 
image  averaging  after  motion  compensation  as  discussed  below. 

The  phase  distortions  suffered  by  individual  pixels  of  an  image  propagated 
through  large  volumes  of  turbulent  atmosphere  may  be  accurately  assumed  zero- 
mean  Gaussian  due  to  the  Central  Limit  theorem.  Furthermore,  the  summation  of 
J  laser-speckle  images,  each  governed  by  an  independent  negative  binomial  process 
with  speckle  parameter  Af,  yields  a  negative  binomially  distributed  composite  image 
with  a  speckle  parameter  of  J  x  Af  [31].  The  composite  image  has  a  mean  intensity 
bound  by  a  much  lower  variance  than  that  of  a  single  speckle  image. 

Image  averaging  reduces  a  large  portion  of  the  random  effects  of  atmospheric  and 
laser  speckle  image  distortion.  However,  the  task  of  image  averaging  is  complicated 
by  atmospheric  tilt  and  platform  motion  or  vibration.  A  hrst  and  necessary  step 
towards  image  averaging  is  multi-frame  registration.  Signihcant  research  has  been 
conducted  towards  the  goal  of  image  registration  in  general,  and  towards  coherent 
image  registration  in  particular  [7,53,62].  The  vector  projection  correlative  regis¬ 
tration  algorithm  of  [7]  is  particularly  attractive  due  to  its  computational  efficiency 
and  accuracy  under  conditions  of  low  SNR.  Since  a  projection-based  algorithm  can 
only  remove  global  tip  and  tilt  in  an  image,  the  resulting  registered  images  will  retain 
distortion  components  from  higher-order  atmospheric  effects  as  well  as  laser  speckle 
degradation.  Figure  2.7(b)  illustrates  the  combined  effects  of  motion  compensation 
and  averaging.  Perfect  image  registration  performance  is  assumed  for  the  analysis  de¬ 
scribed  in  Chapters  III  and  IV,  although  it  is  recognized  that  additional  unmodeled 
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(a)  Simulated  single-frame  image  of  a  tank  (b)  Motion-compensated  frame  averaged  tank 

image 

Figure  2.7:  Subimage  (a)  shows  a  simulated  single-frame  image  of  a  tank  illumi¬ 

nated  by  a  coherent  illuminator  with  a  speckle  parameter  of  M  =  60  and  propagated 
through  a  10  km  turbulent  path  with  a  spherical  tq  value  of  4  cm.  Subimage  (b) 
shows  a  composite  image  created  by  motion  compensating  and  averaging  40  frames  of 
speckle  images  typical  of  those  shown  in  subimage  (a).  Note  the  dramatic  reduction 
of  image  speckle  and  loss  of  high-frequency  image  detail. 


blur  components  will  likely  exist  in  simulated  or  experimental  image  data.  Consid¬ 
eration  is  given  in  Chap.  V  to  the  case  where  image  registration  might  not  provide 
perfect  tip/tilt  compensation  and  thereby  eliminate  motion  blur  due  to  translational 
shift  components  of  the  ensemble  imagery. 

The  average  short  exposure  transfer  function  has  been  shown  to  model  the 
atmospheric  speckle  image  effects  of  a  series  of  tilt-removed  short  exposure  images 
quite  well.  Under  the  assumption  that  individual  speckle  images  collected  by  the 
laser  imaging  system  can  be  regarded  as  independent  realizations  of  images  collected 
through  a  homogeneous  and  isotropic  atmosphere,  the  following  OTF  may  be  used  to 
model  the  motion-compensated  (tip/tilt-removed)  ensemble  average  image  [31], 


Hse  M  =  exp 


(2.39) 


2-31 


where  /  is  the  effective  system  focal  length,  z/  is  a  radial  spatial  frequency  variable,  D 
is  the  diameter  of  the  optical  aperture,  and  ro  is  Fried’s  optical  seeing  parameter.  For 
hxed  optical  components,  the  short  exposure  OTF  is  completely  parameterized  by 
tq.  The  short  exposure  OTF  therefore  provides  an  excellent  candidate  for  use  as  the 
kernel  to  deconvolve  a  blurred  image  in  order  to  obtain  the  original  image  together 
with  the  atmospheric  seeing  condition  parameterized  by  tq. 

For  an  optical  system  that  processes  coherent  illumination,  the  held  propaga¬ 
tion  from  plane  to  plane  may  be  considered  a  linear  process.  However,  the  detected 
intensity  does  not  follow  a  linear  relationship  with  the  originating  intensity  due  to 
the  constructive  and  deconstructive  phase  summations  of  the  random  phasors.  This 
phenomenon  provides  mathematical  insight  towards  the  observation  of  laser  speckle 
in  a  coherent  imaging  system.  In  contrast,  if  a  large  number  of  coherently  formed 
intensity  images  are  averaged  over  a  period  of  time  much  larger  than  the  coherence 
time,  the  composite  image  is  essentially  an  incoherently  formed  intensity  map  of  the 
remote  scene,  and  does  follow  a  linear  relationship  with  the  originating  intensity.  The 
optical  system  may  thus  be  considered  a  linear,  shift-invariant  system  and  the  inten¬ 
sity  of  an  image  pixel,  i{x,y),  may  be  modeled  as  the  convolution  of  the  statistically 
averaged  system  point  spread  function  with  the  remote  scene, 

hysi^^v)  =  X  Hoptiu^v)},  (2.40) 

where  is  the  inverse  Fourier  transform  operator  and  'Hopt{u,v)  is  the  OTF  of 
the  hxed  optical  system.  It  is  assumed  that  the  effects  of  hnite  pixel  size  are  quite 
negligible  compared  to  the  blur  induced  by  the  atmosphere.  With  perfect  motion  com¬ 
pensation,  a  pixel  of  the  average  image  may  be  expressed  as  the  discrete  convolution 
of  hsys{^,v)  with  the  remote  scene,  o{^,ri), 

N  N 

Kx,  y)  =  -i^y-  h)o(C,  y)-  (2.41) 

«=i  »z=i 
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with  images  formed  by  a  detector  of  size  of  iV  x  iV  elements. 


2.9  Knife-Edge  OTF  Estimation  from  Coherent  Imagery 

The  impulse  response  of  a  helded  optical  system  may  be  deduced  from  a  set 
of  experimentally  collected  image  data  from  a  remote  scene  containing  a  step  target. 
Such  a  target  contains  sufficient  contrast  to  extract  the  spatial  step  response  of  the 
system  from  one  of  its  edges.  The  methods  described  in  [78]  were  used  to  recover 
an  estimate  of  Fried’s  seeing  parameter  from  experimental  imagery  collected  by  the 
brass-board  system  described  in  Sec.  2.4. 

A  large  set  of  short-exposure  images  of  a  step  resolution  target  was  collected  by 
the  experimental  laser  vision  system  described  in  Sec.  2.4.  These  data  were  analyzed  to 
produce  an  estimate  of  the  actual  atmospheric  seeing  conditions.  The  individual  image 
frames  were  spatially  registered  by  correlation  with  a  synthetically  generated  step 
target.  This  process  allowed  accurate  motion  compensation  of  the  image  ensemble. 
The  remaining  image  blur  was  then  analyzed  to  estimate  the  seeing  conditions  of  the 
atmosphere  for  the  experimentally  collected  data. 

The  first  step  in  the  process  involves  registering  and  temporally  averaging  hun¬ 
dreds  of  frames  of  image  data  containing  the  step  target.  The  registration  process 
was  executed  on  a  small  portion  of  the  frame  containing  10  pixels  in  the  vertical 
direction  and  20  pixels  in  the  horizontal  direction.  This  configuration  was  used  so 
that  anisoplanatic  effects  would  be  minimized  due  to  the  small  angular  extent  of  the 
target.  Anisoplanatic  effects  cause  increased  blurring  in  the  temporal  average  due  to 
spatially  uncorrelated  motion  which  would  tend  to  bias  the  estimated  seeing  condition 
towards  an  unfairly  low  value  of  tq.  The  non-square  size  of  the  frame  was  chosen  to 
increase  averaging  in  a  dimension  that  did  not  affect  the  registration  process,  since 
the  step  target  contained  no  additional  features  in  the  image  region  where  the  step 
was  not  in  transition. 
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Figure  2.8:  Typical  image  of  the  step  target  used  to  compute  the  impulse  response 
(left)  and  an  image  of  the  ideal  step  target  used  to  register  the  frames  (right). 

As  an  example,  Fig.  2.8  shows  a  typical  frame  of  step  target  data  as  well  as 
the  synthetic  step  used  to  accomplish  the  registration  process.  This  dataset  included 
300  frames  of  imagery  collected  from  a  remote  step-contrast  target  imaged  under 
atmospheric  turbulence  conditions  thought  to  range  from  approximately  1  to  20  cm. 

After  registration  and  averaging,  the  spatial  gradient  in  the  horizontal  direction 
was  computed  in  order  to  estimate  the  derivative  of  the  step  response.  The  derivative 
of  the  step  response  is  the  impulse  response  of  the  system  in  the  horizontal  direc¬ 
tion  [78].  The  short  exposure  impulse  response  was  computed  for  different  values  of 
ro  between  1  and  20  centimeters  in  increments  of  0.1  centimeters  using  the  model 
described  in  Eqn.  2.39  and  a  diffraction  limited  point  spread  function  convolved  with 
a  pixel  of  the  appropriate  size  [60].  Using  an  example  dataset,  it  was  found  that  the 
lowest  mean  squared  error  between  modeled  and  measured  point  spread  functions  was 
obtained  for  a  value  of  tq  =  3.9  cm.  Figure  2.9  shows  the  recovered  impulse  response 
and  the  simulated  impulse  response  for  an  tq  of  3.9  cm.  This  particular  dataset  will 
be  used  in  Chap. Ill  to  provide  atmospheric  truth  to  the  experimentally  collected  data 
used  to  test  the  proposed  blind  deconvolution  algorithm. 
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Figure  2.9:  Impulse  response  recovered  via  knife  edge  calculations  from  the  step 

target  measurements  (solid  line)  and  an  impulse  response  calculated  using  a  short 
exposre  OTF  with  an  Tq  of  3.9  cm  (dashed  line). 


2-35 


III.  Image  Reconstruction  and  Seeing  Condition  Estimation 

using  MAP  Estimation 


The  purpose  of  this  chapter  is  to  develop  a  mathematical  model  of  a  robust  and 
accurate  estimation  algorithm  that  recovers  both  the  remotely  imaged  scene 
intensity,  as  well  as  the  atmospheric  seeing  conditions  under  which  the  imagery  was 
collected.  The  maximum  a  priori  (MAP)  estimator  seeks  a  solution  to  the  blind 
deconvolution  problem  where  an  image  is  presented  to  the  image  processor,  yet  the 
original  truth  scene  and  blur  function  of  the  optical  system  is  not  known.  The  con¬ 
struction  of  such  an  estimator  is  simplihed  somewhat  by  assuming  a  general  form 
of  the  blur  function.  Under  such  an  assumption,  the  research  falls  into  the  more 
restrictive  category  of  parameterized  blind  deconvolution  for  image  reconstruction. 

Active  coherent  illumination  of  remote  scenes  adds  considerable  flexibility  to 
the  task  of  optical  battleheld  sensing.  Unfortunately,  high-power  laser  illuminators 
typically  have  fairly  long  coherence  times,  resulting  in  highly  coherent  scene  illumi¬ 
nation.  Although  multi-frame  averaging  helps  to  mitigate  the  speckled  appearance 
of  the  composite  image  as  discussed  in  Sec.  2.8,  the  statistics  of  the  resulting  image 
cannot  be  accurately  modeled  using  a  Poisson  random  process  unless  relatively  inco¬ 
herent  laser  sources  are  used,  or  alternatively,  a  very  large  number  of  data  frames  are 
combined. 

Images  collected  by  a  partially  coherent  LIDAR  system  experience  atmospheric 
distortion  due  to  the  highly  turbulent  atmosphere  of  the  expected  operating  environ¬ 
ment.  Airborne  imaging  systems  are  subject  to  severe  atmospheric  distortion  due  to 
the  long  slant-range  path  required  of  typical  tactical  scenarios,  while  ground-based 
systems  will  operate  over  shorter  ranges,  but  through  extremely  turbulent  conditions 
caused  by  proximity  to  the  ground.  In  either  case,  typical  seeing  conditions  will  be 
limited  not  by  the  optical  aperture,  but  rather  by  the  atmospheric  coherence  diameter 
quantified  by  Fried’s  seeing  parameter  tq. 

The  development  of  a  joint  estimator  to  discover  both  the  remote  image  re¬ 
flectance  and  atmospheric  seeing  condition  falls  into  the  general  research  category  of 
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coherent  image  reconstruction.  Past  efforts  involving  LIDAR  image  reconstruction 
may  be  roughly  divided  into  adaptive  optic  (AO)  or  post-processing  techniques,  al¬ 
though  there  have  been  some  composite  approaches,  e.g.  [69].  AO  techniques  [60,79] 
promise  good  image  reconstruction  at  the  expense  of  relatively  heavy,  complex  sys¬ 
tems  often  unsuitable  for  manned  and  unmanned  hghter  and  reconnaissance  platform 
deployment.  Furthermore,  such  an  approach  often  requires  point  source  illuminators 
or  guide  stars  to  estimate  the  instantaneous  atmospheric  OTF.  Research  geared  to 
the  problem  of  blind  image  deconvolution  may  be  applied  to  incoherent  image  recon¬ 
struction,  and  such  efforts  are  well  surveyed  in  [41]  and  [42].  However,  the  treatment 
of  blind  deconvolution  of  non- Poisson  distributed  image  sets  is  not  well  covered  in  the 
literature,  nor  has  the  problem  of  parameterized  blind  convolution  been  thoroughly 
studied. 

This  chapter  is  organized  as  follows.  Sections  3.1  and  3.2  describe  the  devel¬ 
opment  of  the  MAP  estimator  for  partially  coherent  multi-frame  image  data,  while 
Sec.  3.2.3  discusses  some  important  implementation  details.  Section  3.3  presents  the 
results  of  estimation  of  seeing  parameters  and  image  recovery  for  both  real  and  simu¬ 
lated  datasets.  Conclusions  and  areas  of  further  research  are  described  in  Section  3.4. 

3.1  Joint  Maximum  a  priori  Image  and  Seeing  Condition  Estimation 

Chapter  II  provided  a  fairly  complete  discussion  of  the  important  details  re¬ 
quired  to  model  a  coherent  optical  system  operating  in  the  presence  of  unknown 
levels  of  turbulence.  The  key  points  of  that  discussion  are  summarized  in  the  con¬ 
text  of  building  an  estimator  that  recovers  an  estimate  of  the  remotely  imaged  scene. 
Without  a  sufficiently  accurate  model  of  the  blurring  effects  due  to  anisoplanatic  tur¬ 
bulence,  the  short  exposure  OTF  provides  a  hrm  foundation  for  parametrization  of 
effects  of  turbulence  on  the  tilt-compensated  optical  system.  The  average  OTF  of 
such  an  optical  system  Ttsys,  including  turbulence,  can  be  mathematically  described 
by 

dtsys  (^7  '^opt  (^5  ("U,  u)  .  (3.1) 


3-2 


where  Tiopt  is  the  non-random  OTF  of  the  optics,  and  Tise  is  the  short  exposure  OTF 
described  by  Eqn.  2.39.  For  hxed  optical  components,  focal  length  and  operating 
wavelength,  the  average  blurring  effects  of  the  atmosphere  are  completely  quantihed 
by  the  atmospheric  coherence  diameter,  or  seeing  parameter,  tq.  This  model  serves 
as  a  reasonable  starting  point  for  the  construction  of  an  image  restoration  algorithm 
under  spatially  invariant  imaging  conditions,  where  the  system  FOV  is  smaller  than 
the  isoplanatic  angle.  Treatment  of  spatially  variant  viewing  conditions  is  deferred 
until  Chap.  IV. 

The  experimental  laservision  system  described  in  Sec.  2.4  captures  a  series  of 
speckle  images  for  post-processing  by  an  off-board  image  restoration  algorithm.  These 
images  are  hrst  spatially  registered  and  then  averaged  to  form  a  motion-compensated 
frame  average  (MCFA)  image  with  reduced  speckle  and  motion  blur.  The  MCFA 
image  is  an  array,  d,  that  represents  a  measurement  of  the  unknown  remote  scene 
o  at  the  imaging  detector.  If  the  characteristics  of  the  atmosphere  at  the  time  of 
image  collection  were  known,  this  information  might  be  used  to  construct  an  estimate 
of  the  short  exposure  OTF  and  could  then  be  used  as  the  deconvolution  kernel  to 
recover  the  best  estimate  of  the  remote  scene  using  an  inverse  Wiener  hlter  or  similar 
techniques  [40].  Unfortunately,  due  to  the  random  nature  of  the  atmosphere  and 
unknown  conditions  likely  to  be  encountered  by  a  fielded  mobile  laser  vision  system, 
an  accurate  estimate  of  the  atmospheric  condition  is  not  usually  available.  Such 
conditions  motivate  the  method  of  blind  deconvolution  by  MAP  estimation  derived 
in  the  following  section. 

3.2  MAP  Estimator  Derivation 

This  section  uses  classical  estimation  theory  to  help  derive  a  joint  MAP  estima¬ 
tor  for  the  remotely  illuminated  scene  together  with  the  atmospheric  seeing  condition 
represented  by  tq. 

Let  D  be  a  random  matrix  representing  a  motion-compensated  frame  averaged 
(MCFA)  image  from  a  collected  ensemble  of  J  speckle  images,  while  O  is  a  non- 
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random  matrix  which  represents  the  remote  scene  or  gronnd  trnth,  Rq  is  a  random 
variable  representing  the  average  spherical  seeing  parameter  and  d,  o  and  tq  are 
specihc  realizations  of  each.  Because  individual  pixel  intensity  has  been  shown  to 
follow  the  distribution  of  Eqn.  2.16,  one  may  express  the  probability  of  the  detected 
image  pixel  given  a  particular  remote  scene  pixel  as 


P[d{x,y)  =  D{x,y)\o{x,y)  =  0{x,y)] 

^  {d{x,y)+M)\  M  ,  i{x,y).M 

{d{x,y)  +  l)\Mr  i{x,y)’  ^  M  ’ 


where  M.  =  J  x  M. frame  is  the  composite  speckle  parameter  of  the  laser  illuminator 
for  the  MCFA  image,  M.  frame  is  the  speckle  parameter  of  each  frame  in  the  ensemble, 
d{x,  y)  is  a  pixel  of  the  MCFA  data,  constructed  from  J  independent  speckle  images, 
and  i{x,y)  is  the  average  intensity  of  a  corresponding  pixel  according  to  Eqn.  2.41. 

Bayes  rule  provides  the  a  posterior  probability  given  the  a  priori  probabilities. 


/o,Ro|d(o,  ro|d) 


/D|o,Ro(d|o,  ro)/RQ(ro)/o(o) 
/D(d) 


(3.3) 


Note  that  the  denominator  of  Eqn.  3.3  is  the  probability  of  a  specihc  realization 
of  a  detected  image.  Although  this  probability  is  not  easily  determined,  it  is  not 
conditioned  upon  the  parameters  of  interest,  and  can  be  treated  as  a  constant  value 
when  forming  a  likelihood  function  [71].  The  prior  /Ro(ro)  is  unknown,  but  may  be 
assumed  as  discussed  below.  The  probability  of  the  object  in  the  numerator,  /o(o), 
is  unknown,  and  may  be  assumed  to  be  a  uniform  distribution.  In  this  case,  Eqn.  3.3 
can  be  more  simply  expressed  as 


/o,iJo|D(o,ro|d) 


/D|o,flo(d|o,  r  o)/flo(r  q) 

/D(d) 


(3.4) 


Although  arbitrary  remote  scenes  certainly  have  some  level  of  spatial  correla¬ 
tion,  image  pixel  distributions  are  assumed  to  be  spatially  independent  from  those  of 
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neighboring  pixels.  The  assumption  of  identically  distributed  and  independent  pixel 
distribution  is  mathematically  convenient  and  common  in  the  derivation  of  maximum 
likelihood  image  estimators  [39,57].  Thus,  the  total  log-likelihood  function  can  be 
expressed  as  the  product  of  the  independent  prior  probabilities 

N  N 

L(o,ro)  =  EEC  ^[fD{x,y)\o(x,y),Ro{d{x,y)\o{x,y),ro)])  ln[//;o(ro)]  (3.5) 

x=l  y=l 

or  by  substitution  with  Eqn.  3.2, 


N  N 


L{o,ro)  =  5^5^  In 


x=l  y=l 


{d{x,y)+M)\ 
[{d{x,y)  +  l)!Af!j 


d{x,  y)  In 


1  + 


M 


Af  In 


1  + 


i{x,y) 

M 


ln[/i?o(A)]- 


(3.6) 


The  probability  density  function,  fn^iro),  represents  the  probability  of  the  seeing 
parameter  tq  being  equal  to  a  specihc  value,  Rq.  The  form  of  the  probability  density 
function  for  the  random  parameter  tq  is  assumed  to  be  exponential  raised  to  the 
power  of  the  number  of  pixels  in  the  array  with  a  mean  determined  by  environmental 
conditions. 


/iJo(A) 


r  _Ar2V2_-| 
e  ’’“”9 


(3.7) 


The  choice  of  this  form  for  the  probability  of  tq  is  based  on  the  empirical  ob¬ 
servation  that  atmospheric  seeing  is  seldom  extremely  better  than  the  average  and 
can  often  be  much  worse.  This  model  also  introduces  numerical  advantages  in  that 
its  logarithm  is  very  simple  to  compute  and  the  entire  distribution  is  characterized 
by  a  single  parameter,  Vavg-  The  influence  of  the  assumed  prior  is  discussed  further 
in  Sec.  3.2.3. 
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3.2.1  Joint  Maximization  of  the  MAP  Likelihood  Funetion.  Joint  parame¬ 
ter  estimation  of  the  remote  scene  together  with  the  seeing  parameter  requires  maxi¬ 
mization  of  an  -|-  1  dimensional  surface.  A  practical  approach  to  maximizing  joint 
likelihood  is  to  select  hxed  candidate  values  for  the  seeing  parameter,  tq  =  Tq,  and  to 
maximize  the  likelihood  of  the  remote  scene  under  the  assumed  atmospheric  condi¬ 
tions.  The  entire  parameter  space  may  be  searched  by  selecting  discrete  values  of  rg. 
Under  this  simplihcation,  the  log  likelihood  of  a  particular  remote  scene  pixel  reduces 
to 

L{o{x,y),rQ)  =  \n[fD(x,y)\o{x,y),Rof{d{x,y)\o{x,y),ro)]  +  \n[fR^{rQ)].  (3.8) 

By  substitution  of  Eqn.  3.8  with  Equations  2.16  and  3.7,  it  can  be  shown  that 


N  N 


{d{x,  y)  \n[i{x,  y)]  -  [d{x,  y)  +  AJ]  ln[i(a;,  y)  +  JV(\) 

x=l  y=l 


-N' 


2  '  0 


—  In 
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avq 

Vjp 


(3.9) 


An  assumption  is  made  that  there  exists  a  particular  realization  of  a  remote 
scene,  6mu  that  maximizes  the  likelihood  of  Eqn.  3.9  at  the  discrete  seeing  condition 
rg,  and  thus,  under  the  necessary  optimality  condition 


(9o(^,r/) 


(3.10) 


Differentiation  of  an  averaged  image  pixel  given  by  Eqn.  2.41  with  respect  to 
the  remote  scene  is  simply 

(3.11) 
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for  each  pixel,  and  hgys  =  T  ^{TLsys)  from  Eqn.  3.1.  The  derivative  of  L(o,  rg)  with 
respect  to  o  is  then 


<9T(o,r[)) 

do{i,ri) 


N  N 

x=l  y=l 


di{x,y) 


{d{x,y)  +M) 


di{x,y) 


M+i{x,y)  do{^,y) 


By  substitution  with  Eqn.  3.11 


<9T(o,r^) 

do{^,v) 


f-f-(d{x,y) 
d{x,  y)  +  M 


hsys{^-x,y-y) 


M  +  i{x,  y) 


hsys{x  -  i,y  -  y)  . 


(3.12) 


(3.13) 


For  the  important  case  where  d{x,y)  =  i{x,y)  for  every  x  and  y,  it  is  easy  to 
show  that  the  functional  is  everywhere  negative,  hence  the  optimality  point  found 
by  Eqn.  3.10  is  a  maximum.  Although  it  is  unlikely  that  direct  solution  to  this 
maximization  problem  might  be  found,  an  iterative  solution  may  be  used.  An  iterative 
maximization  process  may  be  realized  by  setting  Eqn.  3.13  to  zero  and  arranging  terms 
to  yield 


N  N 


,  ,  V  iix,y) 

i;=l  y=l  \ 
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x=l  y=l 


M  +  i{x,y) 


(3.14) 


In  a  manner  consistent  with  the  derivation  of  the  Richardson-Lucy  algorithm 
[57],  both  sides  of  the  equation  may  be  multiplied  by  the  remote  scene  o,  and  an 
update  equation  may  be  formed  to  produce  an  iterative  solution  for  the  estimated 
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scene  parameter 


Z-^x=l  Z-^y=l 


^!^)hsys{x  -  i,y  -  T]) 


EN 

x=l  Z-^y=l 


N  I  d(x,y)+J\4 


with 

N  N 

»”“(x.y)  =  EE  hsysis^  *^5  y 

^-1  ’?=! 


(3.15) 


(3.16) 


In  all  cases  examined,  the  algorithm  produces  a  solution  that  yields  an  estimate 
of  the  observed  scene  given  a  MCFA  image  degraded  by  speckle  and  photon  noise 
at  discrete  values  of  r^.  A  search  over  a  range  of  is  performed  to  hnd  the  MAP 
estimate  of  the  scene  and  Fried’s  atmospheric  seeing  parameter.  As  will  be  discussed 
in  Section  3.3,  the  free  parameter  Vavg  described  in  Eqn.  3.7  seems  to  present  ambi¬ 
guity  in  the  recovered  estimate  of  r^.  It  was  found  that  successive  iterations  of  the 
iterative  algorithm,  starting  with  an  initial  value  of  Vavg  equal  to  the  diameter  of  the 
optical  system  aperture,  provide  an  estimated  seeing  condition  which  was  found  to 
decrease  towards  and  stabilize  on  a  hnal  solution  for  an  estimate  of  r^.  Although 
direct  repetitive  iteration  would  be  time  consuming,  a  more  efficient  approach  will  be 
discussed  towards  the  end  of  Sec.  3.2.3. 

3.2.2  Extension  of  MAP  Estimator  to  Large  Erame  Averages.  It  is  infor¬ 
mative  to  derive  the  MAP  estimator  of  Eqn.  3.15  for  the  case  where  the  number  of 
registered  speckle  images  in  the  ensemble  increases  without  bound.  Alternatively,  the 
coherence  time  of  the  laser  illuminator  may  be  decreased  to  a  point  where  the  speckle 
parameter  becomes  very  large.  In  both  situations,  the  negative  binomial  distribution 
characteristic  of  partially  coherent  illumination  tends  to  the  more  familiar  Poisson 
distribution  that  describes  incoherent  scene  illumination.  Using  a  similar  derivation 
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as  outlined  above,  the  MAP  estimator  becomes 


Z-^x=l  Z-^y=l 


■^^^hsysix  -ty  -v) 


EN  s;^N  7 
X=1  2^17  =  1  ^ 


sys 


(3.17) 


where  the  denominator  acts  as  a  normalization  constant  to  allow  conservation  of  in¬ 
tensity  for  each  iterated  image,  and  i°^^{x,y)  is  given  by  Eqn.  3.16.  This  result  is  an 
expression  of  the  Richardson-Lucy  algorithm,  commonly  applied  to  image  reconstruc¬ 
tion  problems  that  involve  Poisson  noise  processes  [57].  The  contribution  of  Eqn.  3.15 
lies  in  its  ability  to  provide  a  decision-theoretic  MAP  estimate  of  the  remote  scene 
for  cases  where  partially  coherent  illumination  is  unavoidable,  as  in  tactical  scenarios 
where  frame  gathering  time  is  critical,  and  high-powered  laser  illuminators  necessarily 
have  correspondingly  low  speckle  parameters. 


3.2.3  Algorithm  Implementation  and  Choice  for  r^v^.  The  update  algorithm 
for  successive  image  frame  iteration  of  Eqn.  3.15  may  be  easily  implemented  in  a  quasi¬ 
realtime  system  by  recognizing  the  double-summation  of  the  pointwise  divided  images 
as  discrete  convolutions.  Using  fast  Fourier  implementations,  iteration  rates  of  10  —  20 
Hz  are  common  using  modern  desktop  personal  computers  operating  on  256  x  256 
pixel  images.  The  initial  scene  iteration  may  be  started  using  a  uniform  image  matrix 
for  with  a  common  intensity  value  equal  to  the  mean  of  the  input  MCEA  image, 
d.  A  faster  hnal  estimated  solution  is  realized  by  setting  the  equal  to  the  MCEA, 
however,  zeros  contained  in  the  MCEA  data  will  prevent  update  at  the  corresponding 
pixel  location  for  the  hnal  estimated  image,  o”®"'. 

The  mean  intensity  matrix,  i°^'^  is  calculated  iteratively  by  circularly  convolving 
the  system  PSF  with  the  previously  iterated  scene  estimate,  Iterations  may  be 
terminated  when  the  mean-squared  difference  between  the  mean  intensity  1°^*^  and  the 
new  image  convolved  with  the  system  transfer  function  hgys  becomes  less  than 
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the  predicted  variance  of  the  negative  binomial  distribntion  per  Eqn.  5.12,  or  when 


N  N  N  N 


x=l  y=l 


x=l  y=l 


Z-^x=l  Z-^y=l 


d{x,y) 


M 


(3.18) 


Accnrate  algorithm  termination  is  enhanced  by  using  only  bright  pixel  regions  to  form 
estimates  of  image  variance  as  discussed  in  Sec.  2.7. 

In  general,  blind  deconvolution  is  an  ill-posed  problem  [41].  The  assumed  prior 
for  tq  presented  in  Eqn.  3.7  distinguishes  this  MAP  estimator  from  a  maximum  like¬ 
lihood  estimator  by  preventing  the  trivial  solution  of  6  =  d  with  fo  =  cxd.  In  such 
a  case,  the  operator  would  simply  be  presented  with  the  MCEA  image  and  informed 
that  the  atmosphere  caused  no  distortion.  The  assumed  prior  effectively  combats 
selection  of  the  trivial  solution  by  forcing  a  slow  decrease  in  total  likelihood  as  tq  is 
increased.  In  some  situations,  the  average  seeing  condition  may  be  well  quantified. 
In  the  case  of  completely  unspecihed  atmospheric  conditions,  the  introduction  of  the 
prior  introduces  the  free  parameter  Vavg-  To  solve  the  estimation  problem  in  these 
cases,  the  algorithm  may  be  initiated  with  a  larger  than  expected  value  for  ravg,  on  the 
order  of  the  system  entrance  aperture  diameter.  The  MAP  estimate  of  fg  may  then 
be  substituted  for  the  next  estimate  of  Vavg,  and  the  algorithm  repeated  to  iteration 
stopping  criteria  when  fo  ~  Vavg-  This  outer  iteration  does  not  require  that  the  entire 
algorithm  be  run  at  each  iteration.  The  first  two  terms  in  Eqn.  3.9  can  be  computed 
as  a  function  of  tq  in  the  first  iteration  and  saved  as  only  the  logarithm  of  the  prior 
changes  as  a  function  of  ravg-  The  logarithm  of  the  prior  as  a  function  of  ro  and  ravg 
can  be  pre-computed  and  stored  in  a  lookup  table.  This  makes  the  implementation  of 
the  iterative  approach  for  Ending  ravg  as  fast  as  the  implementation  of  the  algorithm 
when  ravg  is  known  a  priori. 


3.3  Results 

This  section  compares  the  MAP  blind  deconvolution  algorithm  using  both  simu¬ 
lated  and  experimentally  collected  data.  Simulated  data  were  constructed  to  compare 
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well  with  resolution  board  and  step  target  data  collected  at  a  mountaintop  test  range 
using  a  candidate  imaging  laser  radar  system  described  in  Sec.  2.4  [54]. 

Table  3.1  describes  the  signihcant  parameters  of  the  simulated  data.  The  exper¬ 
imental  laser  illuminator  speckle  parameter  was  estimated  according  to  the  technique 
documented  in  Sec.  2.5.  Although  12-bit  image  quantization  effects  tended  to  bias 
results  obtained  from  observation  of  dark  pixels,  a  nominal  estimate  of  A4 frame  =  60 
was  obtained  from  observation  of  the  bright  regions  of  several  image  sets.  Motion 
compensated  ensembles  of  J  averaged  images  yielded  composite  speckle  parameter 
estimates  modeled  hj  M.  =  J  x  M. frame  described  in  Sec.  2.8. 

Ground  truth  imagery  was  not  available  in  the  case  of  experimental  image  collec¬ 
tion.  In  this  case,  a  comparison  is  made  between  the  MAP  estimated  seeing  condition 
and  the  seeing  condition  estimated  using  the  knife-edge  line-spread  function  estima¬ 
tion  technique  described  in  Sec.  2.9. 


Table  3.1:  Table  describing  the  simulation  parameters  used  to  create  the  turbulence 
degraded  imagery  used  to  recover  Fried’s  seeing  parameter  using  MAP  estimation. 


Parameter 

Value 

Slant  Range  to  Target 

10  Kilometers 

Optical  Diameter 

20  Centimeters 

Number  of  Phase  Screens 

10 

Distance  Between  Phase  Screens 

1  Kilometer 

Speckle  Parameter  of  Source 

60 

Pixels  per  Image 

128  by  128 

Pixel  Pitch  of  Detector 

11.8  Micrometers 

Mean  Wavelength 

1.54  Micrometers 

Focal  Length 

3  Meters 

Images  per  Frame  Ensemble 

50 

Size  of  Imaged  Target  Area 

5  by  5  Meters 

3.3.1  Results  Obtained  using  Simulated  Image  Data.  Simulated  imagery  was 
created  using  Rayleigh-Sommerfeld  wide  FOV  imaging  as  described  in  Sections  2.1 
through  2.3.  In  order  to  properly  account  for  laser  speckle  effects,  60  images  were 
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Figure  3.1:  Diffraction  limited  simulated  image  reflectance  pattern  used  to  test 

algorithm  performance. 

propagated  through  identically  distributed  but  independent  random  phase  screens. 
The  resulting  images  were  averaged  in  order  to  simulate  a  single  image  with  the  proper 
speckle  parameter  of  Af  =  60.  The  synthetically  generated  image  data  were  quantized 
to  10-bit  resolution  in  order  to  match  the  limited  quantization  of  the  experimentally 
collected  data  described  in  Sec.  3.3.2.  Atmospheric  distortion  was  varied  by  simulating 
various  spherical  tq  values  representing  D/tq  values  of  4,  2.5,  2,  1.3  and  1.0.  One 
thousand  frames  of  partially  coherent  image  data  were  generated  for  each  of  the  five 
atmospheric  conditions.  50-Frame  motion-compensated  frame  average  images  were 
constructed  and  introduced  to  the  iterative  MAP  estimator,  using  uniform  initial  scene 
estimates  with  the  average  seeing  condition  initially  set  to  Vavg  =  D,  the  diameter  of 
the  optics.  The  diffraction-limited  object  is  shown  in  Fig.  3.1. 

Algorithm  iteration  was  usually  complete  within  30-60  seconds  using  a  general 
purpose  PC  running  at  2.8  GHz,  with  differential  image  variance  decreasing  to  the 
analytically  predicted  value  of  Eqn.  2.37  within  approximately  200-300  iterations  for 
low  values  of  rg,  and  approximately  30-50  iterations  as  Tq  values  closer  to  D  were 
searched.  The  change  in  log-likelihood  was  found  to  monotonically  increase  in  all 
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Figure  3.2:  Plot  of  the  number  of  iterative  solutions  required  to  allow  the  estimated 
ro  to  be  within  10%  of  the  hnal  estimated  solution.  Initial  guess  of  Vavg  was  set  to 
the  optical  diameter. 

cases,  with  no  tendency  to  decrease  for  subsequent  iterations  at  a  particular  value 
of  Tq.  Without  introducing  a  stopping  criteria  for  the  iterative  algorithm  at  each 
value  of  Tq,  iterative  estimation  of  the  scene  was  found  to  continue  beyond  the  point 
where  the  recovered  scenes  accurately  represented  the  initial  MCFA  used  as  input  to 
the  algorithm.  Iteration  to  a  stable  value  of  tq  given  an  initial  guess  of  Vavg  =  D 
was  fast  for  scenes  created  with  low  atmospheric  turbulence,  although  the  number  of 
required  iterative  solutions  increased  for  more  turbulent  conditions.  Figure  3.2  shows 
the  number  of  solutions  required  to  move  the  estimated  value  of  fo  from  the  initial 
estimate  of  Vavg  =  D  io  a.  value  within  10%  of  the  hnal  iterated  value. 

Visual  inspection  of  the  resulting  images  indicated  improved  spatial  resolution. 
Table  3.2  shows  the  MAP  estimated  values  of  spherical  %  versus  the  actual  seeing 
parameter  used  to  create  the  turbulence.  Figure  3.4(a)  shows  a  representative  MCFA 
image  for  the  condition  of  H/ro  =  4  (ro  =  5  cm)  used  as  input  to  the  algorithm,  while 
Fig.  3.4(b)  shows  the  recovered  image  using  the  MAP  blind  deconvolution  process. 
Note  that  the  estimated  seeing  conditions  presented  in  Table  3.2  indicate  a  consis¬ 
tently  pessimistic  recovery  of  the  actual  seeing  conditions  used  to  create  the  image 
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data.  It  is  possible  that  such  bias  may  be  attributed  to  motion  blur  effects  induced 
by  anisoplanatic  viewing  conditions.  This  result  is  further  discussed  in  Sec.  3.4. 


Table  3.2:  Simulated  truth  and  estimated  values  for  Fried’s  seeing  parameter  tq 

as  estimated  using  the  MAP  estimation  blind  deconvolution  algorithm  described  in 


Sec.  3.2. 


Simulated  rg  (cm) 

estimated  fp  (cm) 

5 

4.7 

8 

7.5 

10 

9.6 

15 

14.4 

20 

19.1 

3.3.2  Results  Obtained  using  Experimentally  Collected  Image  Data.  Exper¬ 
imentally  collected  data  was  limited  to  a  pair  of  300-image  datasets  collected  for  a 
particular  atmospheric  condition  on  a  controlled  mountain-top  optical  range.  A  reso¬ 
lution  bar  target  and  a  step-intensity  target  were  imaged  according  to  the  parameters 
outlined  in  Table  3.1.  Both  targets  were  arranged  such  that  only  slight  azimuth  change 
was  required  to  image  either  target,  ensuring  similar  atmospheric  prohles.  Addition¬ 
ally,  imaging  of  the  target  sets  was  separated  in  time  by  approximately  2  minutes. 
The  10  km  optical  path  to  the  remote  target  prevented  accurate  atmospheric  truth 
using  scintillometer  measurements.  In  order  to  compare  results,  the  atmospheric  see¬ 
ing  condition  of  the  experimentally  imaged  step-intensity  target  was  estimated  using 
the  knife-edge  OTF  estimation  technique  described  in  Sec.  2.9.  The  averaged  wind 
speed  was  recorded  in  excess  of  35  meters  per  second  at  the  optical  aperture. 

The  experimental  system  used  to  collect  image  data  for  this  effort  was  not  a 
photon-counting  system.  A  calibration  factor  of  p  =  6  photons  per  count  [47]  was 
introduced  to  properly  scale  the  image  intensity  as  discussed  in  Sec.  2.8.  In  addition, 
the  experimental  collection  system  quantized  the  image  intensity  data  into  12-bit 
words.  Dynamic  range  considerations  of  the  analog-to-digital  conversion  process  fur- 
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(a)  Simulated  motion-compensated  frame  aver-  (b)  Recovered  image  using  MAP  estimation 
age  (MCFA)  image  for  blind  deconvolution  algo¬ 
rithm 

Figure  3.3:  Comparison  of  simulated  motion  compensated  frame  average  (MCFA) 
image  under  atmospheric  conditions  of  D/ro  =  4  (ro  =  5  cm)  against  the  MAP 
estimated  image.  Subhgure  (a)  is  the  original  MCFA  image,  while  Subhgure  (b) 
shows  the  MAP  recovered  image  which  was  produced  for  a  most  likely  estimate  of 
ro  =  4.7  cm.  Note  the  additional  high-frequency  spatial  detail  of  Subhgure  (b)  relative 
to  Subhgure  (a). 
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(a)  Experimental  motion-compensated  frame  av¬ 
erage  (MCFA)  image  for  blind  deconvolution  al¬ 
gorithm 


(b)  Recovered  experimental  image  using  MAP 
estimation 


Figure  3.4:  Comparison  of  experimentally  collected  motion  compensated  frame 

average  (MCFA)  to  the  MAP  estimated  image.  Subfigure  (a)  is  the  original  MCFA 
image,  while  Subhgure  (b)  shows  the  MAP  recovered  image  which  was  produced  for  a 
most  likely  estimate  of  tq  =  3.6  cm.  Note  the  additional  high-frequency  spatial  detail 
relative  to  Subhgure  (a).  Estimated  image  contrast  is  somewhat  reduced  by  image 
edge  effects  of  the  deconvolved  image  due  to  the  circular  convolution  properties  of  the 
discrete  Fourier  transform. 


ther  decreased  the  quantization  to  approximately  10  bits  overall  for  the  images  used 
in  this  study. 

Blind  deconvolution  of  both  the  step  intensity  and  resolution  bar  targets  was 
performed,  yielding  an  estimate  of  3.6  cm  for  the  spherical  seeing  parameter,  together 
with  maximum  likelihood  images.  Figure  3.4  (a)  shows  the  10  km  MCFA  image  of 
the  resolution  bar  target  as  input  to  the  algorithm,  while  Fig.  3.4  (b)  shows  the  result 
of  deconvolution.  The  estimate  of  atmospheric  condition  is  in  fairly  good  agreement 
with  a  hgure  of  3.9  cm  derived  from  the  knife-edge  response  of  the  step-intensity 
target. 
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3-4  Conclusions  and  Discussion 

The  agreement  of  estimated  seeing  conditions  with  simulated  and  experimen¬ 
tal  truth  is  encouraging.  Imagery  produced  by  the  automatic  algorithm  consistently 
yields  an  appreciable  increase  in  high  frequency  image  detail,  while  avoiding  the  ten¬ 
dency  to  settle  on  a  trivial  solution  where  the  estimated  scene  equals  the  motion 
compensated  averaged  image  data. 

For  the  case  of  simulated  data,  it  was  noted  that  the  estimates  of  fo  were  approx¬ 
imately  6%  low,  representing  pessimistic  seeing  conditions.  This  discrepancy  might 
be  caused  by  the  relatively  high  levels  of  anisoplanatism  inherent  in  the  simulated 
data. 

Although  the  limited  experimental  data  was  in  close  agreement  with  truth  de¬ 
duced  from  knife-edge  response  techniques,  it  was  noted  that  the  estimate  of  fo  was 
approximately  8%  low.  The  knife-edge  response  technique  was  performed  on  a  small 
region  of  the  step-target  image,  thus  negating  the  additional  blur  caused  by  anisopla¬ 
natism.  Given  conhdence  in  the  truth  data  yielded  by  the  knife-edge  techniques,  it 
appears  that  the  MAP  estimator  is  estimating  poorer  seeing  conditions  than  actually 
encountered  during  data  collection,  since  the  MAP  estimator  uses  the  entire  (aniso- 
planatic)  MCFA  as  input  data.  For  this  particular  viewing  geometry,  it  can  be  shown 
that  the  isoplanatic  angle  is  much  smaller  than  the  held  of  view.  Using  [60],  the 
isoplanatic  angle  is  approximately  24  microradians  for  a  spherical  tq  =  4  cm,  while 
the  system  held  of  view  is  slightly  greater  than  50  milliradians.  However,  the  tilt  iso¬ 
planatic  angle  is  signihcantly  larger  than  24  microradians  [61].  Clearly,  anisoplanatic 
ehects  warrant  consideration  for  these  data. 

The  myriad  of  isoplanatic  patches  that  comprise  the  detected  image  tend  to 
cause  additional  blurring  not  modeled  by  the  short  exposure  OTF  of  Eqn.  2.39,  despite 
ideal  global  image  registration.  A  better  system  model  would  include  the  additional 
blurring  ehects  that  are  produced  by  the  anisoplanatic  imaging  process.  The  analysis 
presented  in  Chap.  IV  attempts  to  resolve  the  quantitative  ehects  of  this  blur,  and 
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replace  the  short  exposure  OTF  with  a  transfer  function  that  captures  the  effects  of 
isoplanatic  patch  motion  blur. 

The  importance  of  accurate  motion  compensation  was  recognized  during  the 
early  portion  of  these  experiments.  Without  accurate  image  frame  registration,  the 
system  models  presented  in  Chap.  II  tend  to  become  less  valid  as  noted  by  inspection 
of  Eqn.  1.1.  In  the  extreme  case  where  no  motion  compensation  is  performed,  the 
application  of  a  long  exposure  OTF  [31]  becomes  appropriate  under  certain  conditions, 
with  the  unavoidable  loss  of  image  detail.  The  problem  of  image  registration  blur  as 
well  as  blur  caused  by  heavy  anisoplanatic  warping  is  considered  in  Chap.  V. 
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IV.  Anisoplanatic  Optical  Transfer  Function  for  Blind 

Deconvolution 


Wide  field-of-view  imaging  systems  present  unique  challenges  to  image  restora¬ 
tion  of  remotely  imaged  scenes.  Common  tactical  employment  of  such  sys¬ 
tems  involves  imaging  paths  though  the  densest  regions  of  the  atmosphere,  over  fairly 
long  horizontal  or  slant  paths  to  target  objects.  In  such  environments,  the  system  FOV 
can  dramatically  exceed  the  tilt  isoplanatic  (isokinetic)  angle,  even  during  modest  lev¬ 
els  of  atmospheric  turbulence  [55].  Under  such  circumstances,  the  statistical  imaging 
model  can  no  longer  be  accurately  described  as  a  shift-invariant  system.  In  contrast, 
the  transmitted  image  is  subject  to  phase  and  amplitude  distortions  that  vary  as  a 
function  of  position  on  the  imaging  device,  due  to  the  wavefront  decorrelation  that 
occurs  due  to  the  necessarily  large  FOVs  required  to  image  typical  scenes-of-interest. 

Several  approaches  have  yielded  successfully  reconstructed  images  under  wide 
FOV  conditions.  A  novel  multiframe  processing  algorithm  is  described  in  [19]  that 
has  been  shown  to  effectively  mitigate  image  degradation  from  coherent  speckle  and 
anisoplanatic  viewing  conditions  by  iteratively  processing  subimage  regions  of  a  re¬ 
mote  scene.  It  appears  that  the  independent  processing  of  multiple  subimages  by  the 
modihed  Ayers-Dainty  blind  deconvolution  algorithm  admits  improvement  for  images 
best  described  by  a  spatially  variant  imaging  process.  Block  matching  or  image  de- 
warping  techniques  [9,23,61],  while  computationally  expensive,  attempt  to  break  the 
shift-invariant  problem  into  several  smaller,  shift-variant  sub-problems.  The  restored 
image  quality  has  been  shown  to  vary  dramatically  as  a  function  of  the  isoplanatic 
patch  size  [9].  However,  without  prior  information  regarding  atmospheric  condition, 
it  becomes  difficult  to  determine  the  number  and  size  of  the  component  patches  that 
comprise  the  detected  image.  Computational  demands  grow  rapidly  as  the  number 
of  isoplanatic  patches  is  increased. 

Adaptive  optic  (AO)  techniques  [33,60,69,79],  may  be  employed  to  mitigate  the 
effects  of  atmospheric  turbulence  in  wide  FOV  systems,  however,  these  techniques  are 
also  computationally  expensive  and  require  substantial  hardware  resources.  In  the 
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case  of  anisoplanatic  viewing  conditions,  multiple  OTFs  corresponding  to  each  iso- 
planatic  patch  must  be  estimated  by  way  of  point-source  illumination  or  artihcial 
guide-star  creation  at  or  near  each  of  the  corresponding  points  on  the  remote  tar¬ 
get  [26,  77].  The  computational  burden  of  multiple  guide-star  creation  and  OTF 
estimation  is  prohibitive  for  application  to  image  reconstruction  from  small,  agile 
tactical  platforms  such  as  manned  hghter-reconnaissance  aircraft  or  remotely-piloted 
vehicles.  While  the  rapid  development  of  micro-electro-mechanical  (MEM)  devices 
will  certainly  revolutionize  the  helding  of  AO  systems,  solutions  that  avoid  AO  archi¬ 
tectures  and  the  associated  computational  burden  are  more  readily  applied  to  space 
and  weight  constrained  applications  in  the  near  term. 

This  chapter  documents  the  derivation  and  application  of  an  anisoplanatic  OTF 
(AOTF)  based  on  tip  and  tilt  correlation  models  of  turbulent  atmosphere  as  the  kernel 
function  of  a  maximum  a  priori  estimation  algorithm  used  to  simultaneously  estimate 
an  image  of  the  remote  scene  together  with  the  atmospheric  seeing  condition  param¬ 
eterized  by  Fried’s  seeing  parameter,  tq.  Previous  research  documented  in  Chap.  Ill 
and  [54]  employed  the  short-exposure  average  OTF,  Hse,  to  model  the  imaging  pro¬ 
cess  of  a  series  of  motion-compensated  speckle  image  frames  from  a  candidate  laser 
vision  system.  Under  this  model,  the  image  formation  process  was  considered  linear 
in  intensity  and  shift-invariant  in  the  average  of  many  such  motion-compensated  im¬ 
age  frames.  However,  it  was  understood  that  T-Lse  was  only  an  approximation  to  the 
true  average  OTF,  since  additional  blur  components  contributed  by  the  motion  of  the 
myriad  of  isoplanatic  patches  were  not  accurately  captured  in  this  statistical  model. 
A  more  accurate  OTF  would  capture  these  additional  blur  components  and  evolve  as 
a  function  of  the  FOV  that  surrounds  the  scene-of-interest. 

This  chapter  is  organized  as  follows.  A  brief  introduction  to  the  difficulties  in¬ 
volved  with  imaging  through  turbulence  with  wide  FOV  systems  is  given  in  Sec.  4.1.1. 
Section  4.1.2  outlines  a  procedure  used  to  create  simulated  anisoplanatic  partially  co¬ 
herent  speckle  imagery  that  is  used  to  compare  to  experimentally  collected  imagery. 
A  simple  tilt-only  model  is  presented  to  describe  the  additional  blur  resulting  from  the 
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random  motion  of  isoplanatic  patches  in  Sec.  4.1.3,  resnlting  in  an  improved  overall 
OTF  to  describe  the  degradation  of  the  motion-compensated  frame  average  image. 
Section  4.1.4  presents  a  model  that  captures  the  additional  tilt  variance  as  a  function 
of  system  FOV  for  a  multi-layered  turbulence  path  to  be  incorporated  into  an  aniso- 
planatic  OTF  for  use  within  a  blind  deconvolution  algorithm.  In  Sec.  4.2,  images  are 
reconstructed  and  seeing  conditions  are  estimated  from  both  simulated  and  experi¬ 
mental  image  ensembles,  using  each  the  original  short  exposure  OTF  and  the  AOTF 
for  comparison.  The  results  are  discussed  and  the  chapter  is  summarized  in  Sec.  4.3. 


4.1  Anisoplanatic  Blur  Model 

4.1.1  Limitations  of  the  Short  Exposure  OTF.  The  FOV  of  tactical  laser 
radar  imaging  systems  is  necessarily  wide,  typically  exceeding  the  isoplanatic  angle  of 
the  atmosphere  by  a  large  margin.  Such  conditions  discourage  modeling  the  formation 
of  images  using  a  single  OTF  due  to  spatial  variance  imposed  by  the  atmosphere. 
However,  it  has  been  demonstrated  that  the  expected  value  of  some  statistical  OTF 
can  accurately  model  the  additional  blur  induced  by  the  uncorrelated  motion  of  the 
multitude  of  isoplanatic  patches  projected  to  the  imaging  detector  [23,69,72].  Such 
an  OTF  must  evolve  as  the  system  FOV  is  changed. 

In  one  extreme,  system  FOV  might  be  made  sufficiently  narrow  that  the  ex¬ 
pected  atmospheric  effects  can  be  accurately  modeled  by  the  average  short  exposure 
transfer  function  introduced  in  Sec.  2.8  and  repeated  here  for  convenience  [31] 


Hse  =  exp 


(4.1) 


where  A  is  the  mean  wavelength  of  the  laser  illumination,  /  is  the  system  focal  length, 
z/  is  a  radial  spatial  frequency  variable,  D  is  the  diameter  of  the  optical  aperture,  and 
ro  is  Fried’s  seeing  parameter  or  optical  coherence  diameter.  For  hxed  optical  compo¬ 
nents,  the  short  exposure  OTF  is  completely  parameterized  by  tq.  The  short  exposure 
OTF  predicts  the  ensemble  average  atmospheric  image  degradation  after  removal  of 
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image  tilt  common  to  the  entire  imaged  scene.  The  OTF  is  accurate  for  an  image 
constructed  from  a  series  of  perfectly  registered  (global  tilt-removed)  image  frames, 
collected  from  a  system  with  such  a  narrow  FOV  that  the  image  transformation  pro¬ 
cess  is  spatially  invariant.  However,  the  short  exposure  OTF  is  optimistic  (predicts 
excessive  high  frequency  detail)  in  the  common  situation  where  system  FOV  begins 
to  exceed  the  isoplanatic  angle  imposed  by  the  turbulent  atmosphere. 

Unfortunately,  tactical  sensors  require  a  very  wide  FOV.  Typical  geometry  con¬ 
straints  of  tactical  sensors  require  that  the  optical  paths  arising  from  individual  points 
that  comprise  an  extended  remote  scene  pass  through  distinct  parts  of  the  turbulent 
atmosphere.  Figure  1.4  depicts  the  geometry  of  a  system  that  experiences  anisopla- 
natic  effects.  Paths  traced  from  a  pair  of  point  sources  separated  by  some  distance  to 
the  telescope  aperture  traverse  regions  of  turbulence  that  possess  different  indices  of 
refraction  and  thus  tend  to  delay  the  optical  phase  by  varying  amounts. 

Excellent  examples  of  relevant  slant-path  propagation  problems  are  presented 
in  [55]  as  well  as  [79]  and  [23],  where  the  effects  of  tilt  anisoplanatism  are  studied 
in  some  detail.  As  a  further  example,  consider  the  following  candidate  laser  vision 
scenario.  An  armored  tank,  with  a  maximum  lateral  extent  of  10  meters  is  viewed 
from  a  distance  of  10  kilometers  using  a  gated  eye-safe  laser  radar  imaging  system. 
The  held  of  view  subtends  approximately  200  arc  seconds  (10  milliradians) .  Further 
assume  a  mean  optical  wavelength  of  1.5  micrometers,  and  a  constant  turbulence  level 
across  an  essentially  horizontal  imaging  path  of  =  10“^^  .  The  isoplanatic 

angle  of  an  arbitrary  optical  system  using  spherical  wave  propagation  is  given  by  [61] 


-3/5 


1 


(4.2) 


where  L  is  the  atmospheric  path  length.  The  resulting  isoplanatic  angle  is  a  mere 
0.55  arc  seconds  (0.28  microradians).  Although  the  tilt  isoplanatic  angle  is  roughly  an 
order  of  magnitude  larger  than  this  hgure  [26] ,  the  system  FOV  remains  dramatically 
larger  than  the  tilt  isoplanatic  angle.  In  wide  FOV  situations  the  number  of  tilt 
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isoplanatic  patches  may  even  approach  or  exceed  the  number  of  pixels  in  the  detector 
array. 


4- 1-2  Simulation  of  Image  Propagation  through  Turbulent  Atmosphere.  Sim¬ 
ulation  of  the  anisoplanatic  imagery  used  for  comparison  to  the  measured  experimen¬ 
tal  data  is  described  in  detail  in  Sec.  2.1  through  Sec.  2.3  but  will  be  briefly  reviewed 
for  clarity.  The  propagation  of  spherical  waves  from  point  sources  that  compose  a 
target  scene  is  disturbed  by  the  random  index  of  refraction  of  turbulent  eddies  within 
the  included  atmosphere  enroute  to  the  imaging  system.  The  large  volume  of  atmo¬ 
sphere  between  the  target  and  optical  aperture  may  be  modeled  as  a  series  of  thin 
phase  screens  placed  at  intervals  along  the  propagation  path.  Due  to  the  large  num¬ 
ber  of  turbulent  eddies  of  varying  refraction  index,  the  central  limit  theorem  permits 
the  phase  delay  of  each  phase  screen  to  be  modeled  as  having  a  circularly  Gaussian 
distribution  of  phase  delay  [35,38].  Phase  screens  are  approximated  as  being  statis¬ 
tically  independent  due  to  physical  separation.  A  single  thin  phase  screen  may  be 
constructed  by  summation  of  the  discrete  individual  screens  while  accounting  for  ge¬ 
ometric  propagation  through  each  screen  [8].  Isoplanatic  effects  of  propagation  may 
be  effectively  simulated  by  conducting  Rayleigh-Sommerfeld  propagation  from  each 
point  on  a  target  through  this  single  thin  phase  screen  placed  at  the  optical  aperture. 
A  coherent  system  model  described  in  Sec.  2.1  and  [8]  is  used  to  create  speckle  images 
at  the  imaging  detector  with  the  correct  spatio-temporal  coherence.  The  composite 
thin-phase  screen  disturbs  independent  realizations  of  speckle  images  in  accordance 
with  the  desired  level  of  atmospheric  turbulence.  Note  that  in  order  to  simulate  aniso¬ 
planatic  viewing  conditions,  a  composite  thin  phase  screen  must  be  created  for  each 
point  propagated  from  the  target  to  the  aperture. 

4-1.3  Optical  Tilt  Effects  Induced  by  Atmospheric  Turbulence.  Tip  and 
tilt  effects  imposed  by  the  turbulent  atmosphere  cause  the  majority  of  image  blur 
in  the  averaged  intensity  of  the  detected  image  ensemble.  Excluding  piston  effects, 
roughly  89%  of  the  turbulence  distortion  power  is  contained  in  Zernike  coefficients  Z2 


4-5 


and  Z3  [60].  The  accurate  spatial  registration  of  a  series  of  short-duration  exposure 
speckle  images  yields  an  average  image  created  by  a  spatially  invariant  system  model 
characterized  by  the  short-exposure  OTF  of  Eqn.  4.1.  However,  in  those  cases  where 
the  atmospheric  viewing  conditions  do  not  permit  system  description  by  a  spatially 
invariant  model,  it  becomes  necessary  to  hnd  other  methods  to  quantify  image  degra¬ 
dation  effects.  A  major  component  of  these  unquantihed  effects  is  due  to  the  blur 
induced  by  averaging  images  which  have  signihcant  amounts  of  local  image  warping 
due  to  anisoplanism.  Another  source  of  blur  in  the  averaged  image  might  be  due 
to  poor  image  registration,  although  this  effect  is  not  treated  in  the  following  analy¬ 
sis  as  it  is  highly  dependent  on  the  performance  of  the  image  registration  algorithm 
employed.  The  focus  of  this  section  is  the  development  of  an  expression  for  an  aniso- 
planatic  OTF  that  captures  the  expected  value  of  the  motion  blur  resulting  from  the 
uncorrelated  tilt  variance  of  point  sources  that  originate  from  the  remote  target  scene 

The  following  analysis  is  presented  to  quantify  the  blur  due  only  to  the  motion 
of  a  point  source  disturbed  by  the  tip  and  tilt  components  of  a  randomly  turbulent 
atmosphere.  Such  analysis  will  allow  an  elegant  description  of  the  blur  that  results 
from  the  decorrelated  motion  between  points  sources  separated  by  greater  than  the 
anisoplanatic  angle.  Considering  the  optical  system  response  to  a  single  target  point 
source  located  close  to  the  optical  boresight,  the  system  may  be  approximated  as  shift- 
invariant.  The  point  spread  function,  p{x,y)  is  dehned  as  the  linear  shift-invariant 
response  of  the  optical  system  to  the  2-D  Dirac  delta  function,  S 

00 

p  {x,  y)=5  (x,  y)  ^  h{x,y)  =  jj  h  (^,  p)  5  {x  -  i,y  -  y)  d^dy  =  h  {x,  y)  (4.3) 

—  00 

where  h  {x,  y)  is  the  instantaneous  impulse  response  of  the  optical  system  and  0 
represents  the  convolution  operator.  To  capture  the  blurring  effects  of  the  atmospheric 
tip  and  tilt  components  exclusively,  one  may  consider  that  the  atmosphere  produces 
a  blur  by  only  shifting  the  location  of  the  point  spread  function  as  a  function  of  time. 
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In  this  case,  the  instantaneous  PSF  is 


p{x,y)  =  5{x  -  a{t)  ,y  -  pit)) ,  (4.4) 

where,  a  and  P  are  random  variables  that  may  be  explicitly  written  as  a  function  of 
time. 

By  definition,  the  Optical  Transfer  Function  is  simply  the  Fourier  transform  of 
the  PSF, 

OO 

T-i  ifx,  fy)  =J^[h  (x,  y)]  =  JJ  p  (x,  y)  (4.5) 

—  OO 

and  by  substitution  with  Eqn.  4.4, 

OO 

ifx,  fy)  =  T[hix,y)]  =  JJ  6  (x  -  ait)  ,y  -  P  it))  e~^^^^^^'^^^^y'>dxdy.  (4.6) 

—  OO 

The  solution  to  the  integral  is  trivial  due  to  the  properties  of  the  Dirac  delta, 

^  ifx,  fy,  t)  =  (4.7) 


Due  to  the  large  volumes  of  distributed  turbulence  between  the  target  and  op¬ 
tical  aperture,  the  tip  and  tilt  or  image  jitter  experienced  at  the  aperture  is  commonly 
assumed  to  be  a  zero-mean  Gaussian  random  process  [25].  Thus,  a  and  P  are  inde¬ 
pendent  Gaussian,  zero-mean  random  variables. 


Pa  (aj  = 


(4.8) 


and 


(/3) 


(4.9) 


The  ensemble  average  OTF  is  easily  found  by  joint  expectation  over  a  and  P, 
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n  (/,,  fy)  =  =  jj  (a,  p)  dad(3.  (4.10) 


A  further  assumption  is  made  that  the  tip  and  tilt  are  mutually  independent,  thus 
their  joint  distribution  is  separable, 


n  (/x,  fy)  =  I I  {a)pf3  {(3)  dad(3.  (4.11) 


By  substitution  of  Equations  4.8  and  4.9  into  Eqn.  4.11  the  OTF  can  be  written 


as 


nu^,f,) 


1  1  ff 

\/2naa  \/2nap  J  J 


g-i27r(/:,o+/j,/3)g 


dad(3, 


(4.12) 


which  is  easily  recognized  as  the  Fourier  transform  of  a  pair  of  jointly  independent 
Gaussian  random  variables.  By  use  of  the  appropriate  Fourier  transform  tables  [31], 
it  can  be  shown  that 


nif.jy) 


^2Tial2Tial 


exp  [-TT  {2TTalf^  +  2nalfl)  ■ 


(4.13) 


Since  the  variances  of  the  tip  and  tilt  and  are  always  positive,  the  average 
OTF  reduces  to 

n  ifx:  fy)  =  exp  [-  {2Ti‘^alfl  +  2Ti‘^alfl)]  ’ 

which  is  a  radially  symmetric  Gaussian  function  due  to  assumed  equal  variance  in  the 
tip  and  tilt  axes.  A  similar  derivation  may  be  found  in  the  literature  [52,80],  and  the 
average  OTF  is  often  expressed  as 


n  (u,  v)  =  exp  [-  (27rV^  (A/u)^  +  27rVJ  (A/u)^)]  ,  (4.15) 
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where  u  and  v  are  spatial  frequency  variables  in  the  aperture  plane,  A/  is  the  optical 
scaling  factor  for  Fresnel  propagation  to  the  detector  located  at  a  focal  length  of  / 
meters  for  a  mean  optical  wavelength  of  A,  and  and  are  tilt  variances  in  the  u 
and  V  directions  respectively. 

H  {u,  v)  can  be  seen  to  have  circular  symmetry  by  letting  u  =  \J {Xfuf'  +  (A/n)^ 
and  by  assuming  equal  variance  power  in  both  the  u  and  v  coordinate  axes.  Under 
these  conditions,  o'q  =  =  <7\  and 

n{iy)  (4.16) 

The  simplihed  OTF  is  completely  parameterized  by  the  axis-combined  tilt  vari¬ 
ance,  a\.  Although  derived  above  for  the  case  of  isoplanatic  imaging,  the  Gaussian 
OTF  provides  the  foundation  for  the  construction  of  an  anisoplanatic  OTF.  If  the 
Gaussian  tilt  variance  that  results  from  an  anisoplanatic  imaging  process  might  be 
derived,  this  variance  may  be  substituted  into  the  above  expression  to  yield  a  suit¬ 
able  AOTF  [72].  The  overall  OTF  employed  to  model  an  ensemble-average  of  well- 
registered  anisoplanatic  imagery  is  simply  the  product  of  the  non-random  OTF  of  the 
optics,  the  OTF  that  embodies  atmospheric  effects  after  global  tilt  motion  compen¬ 
sation  and  some  OTF  that  corresponds  to  the  blur  introduced  by  the  uncorrelated 
motion  of  the  many  isoplanatic  patches, 

'Hsys  {u,  v)  =  Hopt  {u,  v)  Hse  (m,  v)  Haotf  (m,  v)  .  (4.17) 

The  last  two  terms  of  Eqn.  4.17  represent  the  expected  OTF  that  describes  the 
atmospheric  imaging  system  model.  The  true  system  model  is  not  spatially  invariant 
for  the  case  of  anisoplanatic  viewing  conditions,  however,  the  system  may  be  consid¬ 
ered  shift-invariant  in  the  average  with  the  inclusion  of  'Haotf  {u,v)  to  account  for 
blur  effects  due  to  anisoplanatism.  A  comment  on  the  relationship  of  this  model  to  the 
short  and  long-exposure  OTF  is  appropriate.  Assuming  no  image  registration  of  indi- 
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Figure  4.1:  Comparison  of  long-exposure,  short-exposure  and  anisoplanatic  OTFs 
with  D/tq  =  5.  The  long-exposure  OTF  models  an  image  created  by  averaging  a  series 
of  atmospherically  distorted  images  without  the  benefit  of  tip  and  tilt  removal.  The 
higher  spatial  frequencies  evident  in  the  short-exposure  OTF  are  a  direct  result  of  such 
motion  compensation  for  images  collected  with  a  spatially  invariant  optical  system. 
In  systems  that  are  spatially  variant  due  to  anisoplanatic  effects,  high  frequency  detail 
is  lost  due  to  the  summation  of  many  uncorrelated  anisoplanatic  patches  within  each 
image,  resulting  in  an  OTF  that  is  conditioned  on  the  degree  of  anisoplanatism. 

vidual  images  within  the  ensemble,  an  accurate  system  model  would  replace  the  last 
two  terms  of  Eqn.  4.17  with  the  long-exposure  OTF  [31].  Such  an  OTF  would  have 
very  little  high  frequency  content,  due  to  considerable  blur  imposed  by  global  motion 
of  each  ensemble  image.  After  global  image  registration,  dramatic  high  frequency 
spatial  detail  is  gained,  however,  not  as  much  as  if  the  viewing  conditions  permitted 
spatially  invariant  image  formation.  Thus,  the  OTF  most  applicable  to  the  actual 
anisoplanatic  viewing  condition  lies  somewhere  between  the  short  and  long-exposure 
OTF,  as  depicted  in  Fig.  4.1. 


4-1 -4  Tilt  Variance  as  a  Function  of  Geometry.  The  instantaneous  displace¬ 
ment  9?  of  an  imaged  point  at  the  imaging  detector  in  units  of  meters  is  the  integrated 
gradient  of  wavefront  phase,  normalized  by  the  area  of  the  system  aperture  [80] 
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(4.18) 


oo 

A/  JJU(r)  V0  dr 

—  OO 

^  =  - 55 - 

2tt  JJ  U  (r)  dr 

—  OO 


where  U  (r)  is  the  non-zero  held  over  the  aperture  extent,  V0  is  the  wavefront  phase 
gradient,  and  r  is  the  2-D  spatial  variable  in  the  plane  of  the  optical  aperture.  In 
the  simple  case  of  linear  gradient  or  tip/tilt,  the  gradient  can  be  replaced  by  an 
appropriate  slope  multiplied  by  the  independent  variable  in  either  cartesian  coordinate 
axis,  u  or  v. 

Considering  a  series  of  uniformly  spaced  thin  phase  screens  as  depicted  in 
Fig  4.2,  the  u-axis  wavefront  tilt  a  at  the  optical  aperture  in  radians  from  an  arbitrary 
imaged  point  P  {u  +  puiV  +  p„)  on  the  target  displaced  from  the  optical  boresight  may 
be  synthesized  from  a  normalized,  weighted  sum  of  basis  vectors 

oo  N 

n  Y.i.u  +  +  pvr,)K{u  +  p^^,V  +  Pvjcpn  {u,  v)  dudv 

-oo  n=l 

«  =  - 55 - 

f  f  uA  (u,  v)  dudv 

—  OO 

where  u  and  v  are  spatial  variables  in  the  aperture  in  units  of  meters  that  correspond 
to  the  X  and  y  coordinates  in  the  detector,  the  random  held  <pn  {u,  v)  represents  the 
phase  of  the  of  N  phase  screens,  and  An{u  +  Pu^^v  +  is  the  deterministic 
aperture  weighting  function  geometrically  formed  by  projection  from  the  point  source 
P  {u  +  pu,v  +  p^)  on  the  target  to  the  optical  aperture  plane,  pu  and  are  the 
orthogonal  components  of  displacement  in  meters  from  optical  boresight  measured  at 
the  target  plane.  Since  the  discrete  atmospheric  phase  screens  are  uniformly  spaced, 
the  displacement  of  the  projected  apertures  from  a  point  Pi{ui,  Ui)  at  each  screen  may 
be  found  according  to  the  linear  relationships  =  ^ui  and  =  ^vi  where  Zn  is 
the  z-axis  position  of  the  phase  screen,  and  L  is  the  distance  from  the  target  to 
the  optical  aperture.  Typical  aperture  weighting  functions  are  radially  symmetrical 


(4.19) 
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Figure  4.2:  Coherent  imaging  model.  Two  distinct  point  sources  from  the  target 

propagate  to  the  optical  aperture  through  a  turbulent  atmosphere  represented  by  N 
thin  phase  screens  uniformly  distributed  along  the  path  from  target  to  receiver.  At 
each  random  atmospheric  thin  phase  screen,  a  projected  aperture  function  is  formed 
by  the  physical  geometry  of  the  point  source  and  receiver  location  as  depicted.  The 
model  is  used  to  predict  the  summed  tilt  contribution  from  each  point  source  at  the 
optical  aperture. 

and  have  a  value  of  unity  within  the  aperture,  and  zero  outside  this  region.  The 
denominator  of  Eqn.  4.19  serves  to  normalize  tilt  magnitude  with  reference  to  the 
aperture  weighting  function  of  the  optical  system,  A  {u,v). 

Heuristically,  the  tilt  measured  at  the  optical  aperture  for  a  given  point  on  the 
target  may  be  understood  to  be  composed  of  two  components.  One  component  is 
a  global  tilt  which  may  be  effectively  removed  by  accurate  image  registration.  The 
second  component  is  local  tilt  due  to  anisoplanatic  viewing  conditions.  Note  that  these 
two  components  are  independent,  since  any  global  tilt  that  might  exist  in  the  local  tilt 
component  would  be  removed  by  global  image  registration.  Therefore,  the  composite 
tilt  a  may  be  written  as  a  =  ac  +  du,  where  dc  and  du  are  the  spatially  correlated 
and  uncorrelated  tilt  components  respectively.  Considering  a  pair  of  arbitrary  points 
on  a  target,  it  may  easily  be  shown  that  their  uncorrelated  tilt  variance  is  the 

residual  variance  calculated  by  subtracting  the  correlated  tilt  variance  from  the  total 
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tilt  variance, 


Ela^]  =  E[a‘^]  —  E[aia2],  (4.20) 

where  i?[Q;^]  is  the  overall  tilt  variance  of  either  of  the  points  at  the  aperture  and 
£^[0:102]  is  the  correlated  tilt  power  between  the  two  points. 

It  is  instructive  to  hrst  compute  the  aperture  tilt  covariance  i?  [0:10:2]  resulting 
from  two  distinct  points  on  the  target,  Pi  {ui,  ui)  and  P2  {u2,  V2).  £'[0102]  will  depend 
only  on  point  separation  if  the  atmospheric  turbulence  can  be  considered  isotropic. 
For  a  fixed  level  of  atmospheric  turbulence,  more  distant  points  on  the  target  will 
produce  less  correlated  tilt  power  at  the  aperture.  On  the  other  hand,  two  closely 
spaced  points,  subtending  an  angle  well  within  the  tilt  isoplanatic  angle  for  a  given 
turbulence  level,  will  yield  essentially  no  uncorrelated  tilt  variance  E[a^. 

Calculation  of  the  tilt  covariance  of  two  arbitrary  points  is  mathematically 
straightforward.  It  is  convenient  to  fix  one  point  at  the  optical  boresight  of  the 
target,  while  locating  the  second  point  at  a  distance  of  pu  and  in  the  u  and  v 
directions  respectively.  The  u-axis  tilt  covariance  at  the  aperture  is  then 


E  [0:10:2] 


tjj 


N  N 

EE 

721=1  n2  =  l 


E 


uA. 


ni 


[U,V] 


[U,V] 


(«'  +  Pnn  J  An2  [u'  +  Pu„^  ,  v'  + 
0n2  ip'  ,v')  dudvdudv'  , 


(4.21) 


where 


^/j  = 


uA  (m,  v)  dudv 


(4.22) 


is  the  scalar  normalization  constant  formed  by  integration  over  the  extent  of  the  u-axis 
tilted  aperture  in  the  plane  of  the  optical  receiver. 
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As  is  common  in  the  literature,  the  individual  phase  screens  are  taken  to  be 
statistically  independent  zero-mean  random  Gaussian  helds,  therefore, 


E  (ui,  ui)  (“2,  ^'2)]  =0,  V  ni  7^  n2 


(4.23) 


which  results  in  the  cancelation  of  summation  cross-terms.  After  dropping  the  sub¬ 
script  on  n  for  convenience,  the  tilt  covariance  becomes 


E  [aias]  /  /  /  /  '^uAn  (u,  v)  {u  +  An  {u  +  p„„,  u'  +  p„„) 


E  [0„  (u,  v)  (pn  u')]  dudvdu'dv' . 


(4.24) 


After  making  the  substitutions  of  variables,  it  =  u'  +  pu„  —u  and  v  =  v'  +  —v, 

the  correlated  component  of  tilt  at  the  aperture  may  be  expressed  as 


E  [0:10:2]  =  /  /  /  /  (m,  v)  [it  +  u)  An{u  +  U,v  +  v) 


E  \4>n  (u,  v)<pn(u+[u-  PuJ,v+  [!'  -  p„,J )]  dudvdudv, 

(4.25) 

and  the  expected  value  operation  in  the  integrand  may  be  recognized  as  an  autocorre¬ 
lation  of  the  phase  screen,  which  is  only  a  function  of  u,v,pu„  and  p„„.  This 
observation  allows  the  quadruple  integral  to  be  split  into  a  pair  of  double  integrals 


E  [01O2]  =  Ip  Eh 


[u-  pun^v-  pv 


uAn  (u,  v)  {it  +  u)  An{u  +  u^v  +  v)  dudv  >  dudv.  (4.26) 
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The  integral  within  the  braces  of  Eqn.  4.26  is  the  2-D  autocorrelation  of  the 
tilted  aperture  function.  Let  this  quantity  be  represented  by  Gu„{u,v).  Then, 


N 


n=l 


E[aia2]=i^Yl  //  ^<Pn{u- Pu„,v  -  PvJGu„{u,v)dudv.  (4.27) 


It  is  now  necessary  to  quantify  the  tilt  variance  of  an  arbitrary  point  in  the 
scene,  E[a‘^].  Similar  analysis  is  presented  in  the  literature,  e.g.  [3].  For  the  case  of 
tilt  variance  of  a  given  point  at  the  target,  there  exists  no  displacement  between  point 
sources  at  the  target,  i.e.,  pu„  =  0  and  =0,  V  n  G  iV,  and  the  uncorrelated  tilt 
at  the  aperture  may  be  found  by  substitution  of  i7[a^]  into  Eqn.  4.20, 


E[al  = 


N  ^ 

V’E  // 

n=l 


{u-0,v  -0)  -  {u-  Pu^,v-  Pyj]  Gu„  {u,  v)  dudv.  (4.28) 


The  structure  function  may  be  expressed  as  (p)  =  (0)  —  (p)) 

for  an  isotropic  turbulent  atmosphere,  where  p  represents  radial  separation  between 
points  on  a  thin  phase  screen.  After  adding  and  subtracting  R^^  (0,  0)  to  both  terms 
within  the  integrand,  the  uncorrelated  tilt  may  be  written 


^  K]  =  ^  E  II  pu„,v  -  pyj  Gu^  {u,  v)  dudv. 


N 


n=l 


(4.29) 


where 


{u-pu^,v-  pyj  =  {u  -  pu^,v  -  pyj}  -  {u,  h)}  .  (4.30) 

An  identical  derivation  can  be  used  to  hud  the  transverse  component  of  tilt. 
A  more  useful  form  of  the  uncorrelated  tilt  power  can  be  found  by  summing  the 
orthogonal  components  of  tilt  variance  in  each  axis.  Assuming  symmetry  of  the 
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optical  aperture  in  both  the  u  and  v  axes,  the  combined  total  uncorrelated  tilt  power 
is  then 


N 


n=l 


^  ^  X]  /  /  {u-  Pu„,v-  p^J  {Gu^  {u,  v)  +  {u,  h)}  dudv.  (4.31) 


By  allowing  pn  =  a/ pI^  +  pl^  and  recognizing  the  radial  symmetry  of  the 
summed  tilted  aperture  correlation  functions,  Gn{u,v)  =  Gu^  {u,v)  +  {u,v),  the 

uncorrelated  tilt  may  be  expressed  as  in  integral  over  polar  coordinates  r  and  6, 


N 

n=l 


Pn)Gn{r)  rdrdO. 


(4.32) 


The  phase  structure  function  noted  in  Eqn.  4.30  may  be  most  simply  modeled 
by  using  Kolmogorov  statistics,  with  the  important  limitation  that  the  effects  of 
turbulence  are  constrained  to  some  inertial  subrange  such  that  turbulent  eddy  sizes 
are  bounded  by  the  upper  and  lower  scale  values  Lq  and  /q  respectively  [1].  Then, 

D4,  (p,  Zn)  =  J  Gl  (z)  l—J  dz,  /o  <  Po  <  Lq.  (4.33) 

Assuming  constant  G^  prohle  as  might  be  encountered  during  a  horizontal  path 
imaging  scenario,  Eqn.  4.33  reduces  to 


(p,  Zn)  =  l.^9Glznk‘^ p^^^ ,  /o  <  Po  <  ^0,  (4.34) 

where  the  spherical  coherence  radius  is  po  =  {Q.bbG‘^Znk‘^)  and  k  =  For  each 
phase  screen,  the  length  Zn  is  the  extent  of  the  atmosphere  encompassed  by  the 
screen  in  the  z-axis. 
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Substituting  the  Kolmogorov  phase  structure  function  from  Eqn.  4.34  into 
Eqn.  4.32  and  expressing  Fried’s  parameter  for  spherical  wave  propagation  as 


r 


O 


47r^ 

^2/^2  3 


(4.35) 


a  form  of  Eqn.  4.32  more  suitable  for  use  within  the  parameterized  blind  deconvolution 
algorithm  is 


N 


E  =  11.6277r^'^r, 


^TT  poo  _ 


0  /  / 


~  {r)  rdrdO,  (4.36) 


where  E  [r^]  is  the  uncorrelated  angular  tilt  variance  at  the  optical  aperture  expressed 
in  units  of  square  radians. 

One  hnal  simplihcation  may  be  realized  due  to  the  radial  symmetry  of  the 
integrand,  with  the  expression  reduced  to  a  single  integral  after  integration  over  P, 


E  [tI]  =  23.2547r^^7r, 


N 


_i  JO 


{r  -  Gnir)dr, 


(4.37) 


n=l 


As  expected,  for  a  target  point  P2  {u2,V2)  located  along  optical  boresight  (i.e., 
\pn\  =  0),  the  uncorrelated  tilt  variance  is  zero,  and  all  of  the  tilt  power  between  the 
two  points  is  correlated.  The  numerically  computed  uncorrelated  tilt  of  a  point  as  it 
is  displaced  from  the  optical  boresight  for  a  particular  anisoplanatic  viewing  geometry 
and  range  of  atmospheric  conditions  is  calculated  according  to  Eqn.  4.37  and  plotted 
in  Fig.  4.3. 

The  radially  symmetric  integrand  of  Eqn.  4.37  that  results  from  the  multiplica¬ 
tion  of  the  auto-correlation  of  the  tilted  aperture  functions  Gn  (r),  with  the  displaced 
structure  function  (r  —  p„)  is  only  a  function  of  the  imaging  system  geometry 
given  a  particular  thin  phase  screens  location,  and  may  be  pre-computed  in  an  off¬ 
line  system  for  various  target  engagement  ranges.  A  fast  look-up  table  approach 
would  obviate  the  need  for  calculation  of  the  2-D  cross  correlations  of  the  tilted  aper- 
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Figure  4.3:  Uncorrelated  tilt  variance  as  a  function  of  angular  point  separation. 

The  total  uncorrelated  tilt  variance  is  plotted  as  a  function  of  point  separation  in 
the  target  plane  for  several  values  of  tq.  The  optical  aperture  is  20  cm,  and  range  to 
target  is  10  km. 

ture  functions  in  a  real-time  environment.  The  numerical  calculation  of  the  aperture 
auto-correlation  Gn  (r)  is  required  only  at  a  single  phase  screen,  and  may  be  linearly 
scaled  for  the  remaining  screens.  As  an  alternative  to  numerical  calculation  of  the 
integral  in  Eqn.  4.37,  (r)  should  be  easily  computed  in  closed  analytical  form, 

given  relatively  simple  aperture  geometries  [63].  Figure  4.4  shows  an  example  tilted 
correlation  function  for  the  common  case  of  a  uniform  circular  aperture  weighting 
function. 

The  well-behaved  nature  of  the  structure  function  in  the  limit  as  |/7„|  approaches 
zero  allows  the  use  of  any  of  the  common  atmospheric  models,  including  Kolmogorov, 
von-Karman,  modihed  von-Karman  [38],  etc.  However,  an  important  limitation  of 
the  expression  when  used  for  turbulence  strength  estimation  is  the  requirement  to 
assume  constant  as  a  function  of  distance  to  the  target.  Such  an  assumption 
allows  compact  parameterization  of  the  seeing  condition  using  spherical  tq  and  is  often 
made  for  horizontal  and  moderate  slant-path  imaging  as  is  typically  encountered  in 
the  tactical  observation  environment. 
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Figure  4.4:  2-D  autocorrelation  of  a  tilted  circular  aperture  weighting  function. 

This  Witch ’s  Hat  function  was  created  by  the  sum  of  the  u-axis  tilted  autocorrelation 
function  with  that  of  the  v-axis  tilted  function. 

y^.i.5  Anisoplanatic  OTF  for  Wide  FOV  Systems.  The  analysis  conducted 
in  Sec.  4.1.4  describes  the  uncorrelated  tilt  power  at  the  optical  aperture  due  to  a  single 
point  source  on  the  target  separated  by  some  distance  \pn\  from  the  optical  boresight. 
Clearly,  larger  FOV  systems  will  be  best  described  by  AOTFs  that  incorporate  greater 
uncorrelated  tilt  variance.  To  achieve  an  appropriate  OTF,  a  snitable  point  separation 
must  be  chosen  to  captnre  the  expected  uncorrelated  blur  effects  that  span  the  system 
FOV.  An  argnment  can  be  made  to  select  a  radial  point  separation  equal  to  the  radius 
of  the  system  FOV.  Such  a  choice  yields  an  OTF  that  predicts  the  maximum  amount 
of  nncorrelated  tilt  motion  blur  in  each  averaged  data  frame, 

hiAOTF  (4.38) 

where 

^  poo 

^  =  23.2547rW^r-5/3  ^  ^  _  ^5/3  (4  39) 

is  obtained  from  Eqn.  4.37  with  pn  =  ppoVr^-  Here,  is  the  tilt  variance  expressed 

in  nnits  of  sqnare  meters  at  the  detector  plane  corresponding  to  the  conical  held  of 
view  described  by  a  point  chosen  at  the  furthest  extent  of  the  target  image  under 
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observation.  The  scale  factor  of  is  required  to  transform  aperture  wavefront  tilt 
variance  to  image  displacement  variance  at  the  detector  plane. 

For  a  hxed  system  focal  length,  an  operational  user  might  choose  to  reduce  Pfov 
in  order  to  surround  only  the  “objects  of  interest”  in  the  overall  target  scene,  allowing 
a  more  appropriate  OTF  for  use  in  a  suitable  image  restoration  algorithm. 

4-2  Results 

A  maximum  a  posteriori  blind  deconvolution  algorithm  [54]  that  estimates  the 
remote  scene  image  together  with  the  atmospheric  seeing  condition  was  used  to  process 
both  simulated  and  experimentally  collected  data.  The  MAP  algorithm  is  briefly 
reviewed  in  Sec.  4.2.2.  As  noted  in  Chap.  Ill,  the  short  exposure  OTF  was  employed 
to  form  a  MAP  estimate  for  Fried’s  seeing  parameter,  tq.  Although  image  detail  was 
restored  relative  to  the  motion-compensated  frame  average  image  data  input  to  the 
routine,  it  was  understood  that  the  additional  blur  introduced  by  anisoplanatic  patch 
motion  was  not  properly  modeled  via  the  short  exposure  OTF  of  Eqn.  4.1. 

This  section  compares  the  recovered  seeing  parameters  and  images  obtained 
using  the  improved  AOTF  to  those  estimated  using  the  short-exposure  OTF.  Sec¬ 
tion  4.2.1  documents  the  procedure  used  to  create  the  simulated  data  and  the  method 
used  to  extract  atmospheric  truth  from  the  experimentally  collected  data.  Sec¬ 
tion  4.2.3  outlines  the  compared  results  using  simulated  wide  FOV  speckle  image 
data,  while  Sec.  4.2.4  documents  the  comparison  for  experimentally  collected  data 
obtained  from  an  optical  range  located  at  North  Oscura  Peak,  White  Sands  Missile 
Range  (WSMR),  New  Mexico. 

4-2.1  Experimental  Method.  To  quantify  the  improvement  of  tq  estimation 
over  a  broad  range  of  seeing  conditions,  simulated  anisoplanatic  10-Kilometer  laser 
radar  data  were  generated  and  processed  using  both  the  short  exposure  OTF  and  the 
improved  anisoplanatic  OTF.  The  simulated  image  data  consisted  of  a  subsection  of 
a  standard  resolution  board  target  with  a  system  FOV  that  greatly  exceeded  the  tilt 
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isoplanatic  angle  of  the  simulated  turbulent  atmosphere  for  various  levels  of  Fried’s 
seeing  parameter.  Experimentally  collected  data  consisted  of  hve  large  ensembles  of 
speckle  imagery  of  a  standard  resolution  target  and  step-contrast  target  located  at  a 
range  of  10  km  from  the  laser  vision  system. 

4. 2. 1.1  Synthetic  Imagery  Generation.  Prior  to  consideration  of  atmo¬ 
spheric  effects,  simulated  image  data  were  generated  to  capture  the  effects  of  partially 
coherent  laser  illumination  of  a  remote  scene.  The  spatio-temporal  coherence  proper¬ 
ties  of  a  gated  laser  illuminator  may  be  effectively  characterized  by  a  scalar  speckle 
parameter,  M..  The  parameter  may  be  mathematically  regarded  as  the  degree-of- 
freedom  parameter  of  the  negative  binomial  distribution  used  to  describe  the  detec¬ 
tion  of  coherent  illumination  at  the  detector  [31].  Note  that  in  the  limit  as  M.  grows 
without  bound,  the  negative  binomial  distribution  tends  toward  the  familiar  Pois¬ 
son  distribution  often  used  to  model  incoherent  illumination.  Physically,  M.  may  be 
understood  to  represent  the  degree  of  illuminator  stability  over  an  observed  gating 
period.  Laser  illuminators  possessing  large  speckle  parameters  have  less  coherence 
and  detected  images  tend  to  have  less  laser-speckle. 

A  perfectly  coherent  plane  wave  was  assumed  incident  on  the  target.  The  target 
coordinates  were  assigned  such  that  the  target  depth  along  the  axis  of  propagation 
varied  uniformly  with  a  variance  of  10  optical  wavelengths.  Such  target  roughness 
caused  optical  phase  interference  and  yielded  a  completely  developed  intensity  speckle 
pattern  at  the  detector  subject  to  the  diameter  of  the  optical  aperture.  Images  pro¬ 
duced  from  a  system  with  shorter  coherence  times  were  simulated  by  averaging  several 
fully  developed  speckle  images.  As  an  example,  a  simulated  image  collected  from  a 
laser  system  found  to  have  a  speckle  parameter  of  Ad  =  60  requires  the  generation 
and  summation  of  60  fully  developed  speckle  images,  each  created  with  different  uni¬ 
formly  distributed  random  surface  roughness.  Table  4.1  lists  the  salient  parameters 
used  to  conduct  the  simulation. 
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Table  4.1:  This  table  describes  the  simulation  parameters  used  to  create  the  aniso- 
planatic  turbulence  degraded  imagery  used  to  recover  Fried’s  seeing  parameter  by  way 
of  blind  image  deconvolution. 


Parameter 

Value 

Slant  Range  to  Target 

10  Kilometers 

Optical  Diameter 

20  Centimeters 

Number  of  Phase  Screens 

10 

Distance  Between  Phase  Screens 

1  Kilometer 

Speckle  Parameter  of  Source 

60 

Pixels  per  Image 

128  by  128 

Pixel  Pitch  of  Detector 

11.8  micrometers 

Mean  Wavelength 

1.54  micrometers 

Focal  Length 

3  Meters 

Images  per  Frame  Ensemble 

50 

Size  of  Imaged  Target  Area 

5  by  5  Meters 

To  capture  the  deleterious  effects  of  the  turbulent  atmosphere,  the  10-kilometer 
propagation  path  was  divided  into  ten,  1-kilometer  atmospheric  volumes.  Statis¬ 
tically  independent  random  Gaussian  thin  phase  screens  were  constructed  using  a 
method  similar  to  that  documented  in  [35]  and  [38]  for  each  of  the  ten  volumes.  The 
phase  delay  effects  for  each  of  the  atmospheric  volumes  were  effectively  collapsed  to 
a  thin  phase  screen  located  behind  each  volume.  To  capture  anisoplanatic  effects, 
the  projected  sub-apertures  for  each  imaged  point  were  calculated  at  each  of  the  ten 
thin  phase  screens.  The  optical  phase  delay  through  each  of  these  sub-apertures  was 
summed  to  create  a  single  thin  phase  screen  located  at  the  optical  aperture.  For  each 
imaged  point  at  the  remote  target,  light  was  propagated  using  Rayleigh-Sommerfeld 
propagation  to  the  optical  aperture  [8],  [54].  To  simulate  the  effect  of  partially  coher¬ 
ent  illumination,  M  fully  developed  speckle  images  were  propagated  using  varying 
surface  roughness  for  each  set  of  phase  screens  and  averaged  to  create  a  single  aniso- 
planatically  distorted,  partially  coherent  speckle  image. 
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4-2. 1.2  Knife-Edge  OTF  Estimation  from  Experimentally  Collected  Im¬ 
agery.  The  following  section  provides  a  review  of  the  knife-edge  techniques  intro¬ 
duced  in  Sec.  2.9.  A  large  image  set  of  a  step  resolution  target  ground  truth  was 
collected  by  the  experimental  laser  vision  system  described  in  Sections  2.4  and  4.2.4. 
These  data  were  analyzed  to  produce  an  estimate  of  the  actual  atmospheric  seeing 
conditions.  The  individual  image  frames  were  spatially  registered  by  correlation  with 
a  synthetically  generated  step  target.  This  process  allowed  accurate  motion  com¬ 
pensation  of  the  image  ensemble.  The  remaining  image  blur  was  then  analyzed  to 
estimate  the  seeing  conditions  of  the  atmosphere  for  the  experimentally  collected  data. 

The  long  path  between  the  imaging  system  and  target  makes  atmospheric  seeing 
condition  measurement  difficult  using  standard  scintillometry  techniques.  To  obtain 
accurate  atmospheric  truth,  a  line-spread  function  was  deduced  from  experimentally 
collected  imagery  according  to  the  method  described  in  [78].  Five  large  sets  of  im¬ 
age  data  were  considered,  each  consisting  of  300  speckle  images.  For  each  data  set, 
the  speckle  images  of  a  step-contrast  target  was  hrst  registered  and  then  averaged 
to  produce  a  single,  motion-compensated  image  frame.  The  spatial  gradient  in  the 
horizontal  direction  was  computed  from  this  image  in  order  to  estimate  the  derivative 
of  the  step  response.  The  derivative  of  the  step  response  is  the  impulse  response  of 
the  system  in  the  horizontal  direction  [78].  The  short  exposure  impulse  response  was 
computed  for  different  values  of  ro  between  1  and  20  centimeters  in  increments  of  0.1 
centimeters  using  the  model  described  in  Eqn.  2.39  and  a  diffraction  limited  optical 
transfer  function  convolved  with  a  pixel  of  the  appropriate  size  [60].  The  simulated 
OTFs  were  best  £t  to  the  data  using  the  least-square  error  metric.  The  estimated 
spherical  seeing  conditions  for  each  data  set  are  tabulated  in  Table  4.2. 

4.2.2  MAP  Blind  Deconvolution  Algorithm.  A  novel  blind  deconvolution  al¬ 
gorithm  was  previously  developed  to  simultaneously  estimate  a  remote  scene  together 
with  current  seeing  conditions  [54].  The  algorithm  was  developed  using  a  Bayesian 
approach  under  the  assumption  that  the  detection  of  partially  coherent  illumina- 
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tion  follows  a  negative  binomial  statistical  distribution  [31].  Under  this  framework, 
a  likelihood-based  cost  function  was  constructed  for  each  pixel  of  the  detected  im¬ 
age.  An  assumption  of  independent  pixel  distributions  was  made,  whereby  the  total 
likelihood  function  for  the  image  was  formed  from  the  product  of  the  individual  dis¬ 
tributions.  The  likelihood  cost  function  was  modihed  by  the  addition  of  an  assumed 
prior  for  the  seeing  condition.  The  assumed  prior  followed  a  negative  exponential 
distribution,  due  to  the  physical  observation  that  atmospheric  seeing  is  seldom  ex¬ 
tremely  better  than  some  average  condition  and  can  often  be  worse.  For  each  value 
of  the  seeing  condition  characterized  by  Fried’s  seeing  parameter,  an  iterative  maxi¬ 
mization  of  the  likelihood  was  performed.  For  each  specihc  value  of  tq,  iterations  were 
continued  until  the  variance  of  the  image  decreased  to  that  predicted  by  the  negative 
binomial  distribution.  It  was  found  that  the  algorithm  revealed  a  maximum  value  of 
ro,  beyond  which  the  likelihood  tended  to  decrease  due  to  the  influence  of  the  negative 
exponential  prior.  This  estimate  for  tq,  together  with  the  resulting  estimated  image, 
represented  a  useful  solution  to  the  blind  deconvolution  problem.  However,  it  was 
understood  that  the  employed  deconvolution  kernel  did  not  account  for  the  effects 
of  atmospheric  anisoplanatism,  and  tended  to  provide  pessimistic  estimates  of  seeing 
conditions. 

4-2.3  Results  Obtained  using  Simulated  Image  Data.  Five,  1000-image  data 
sets  were  constructed  to  simulate  atmospheric  conditions  described  by  spherical 
values  of  5,  8,  10,  15  and  20  centimeters,  for  D/vo  values  of  5,  2.5,  2,  1.33  and 
1.0  respectively.  A  MAP  blind  deconvolution  algorithm  [54]  was  used  to  estimate  the 
most  probable  seeing  parameter  for  each  of  the  averaged  images  formed  from  50-frame 
ensembles  within  each  data  set.  The  algorithm  was  run  using  both  the  short-exposure 
OTF  and  the  new  AOTF.  Iteration  was  allowed  to  continue  until  the  variance  between 
the  convolved  estimated  image  decreased  to  the  variance  of  the  negative  binomial 
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distribution  used  to  model  the  statistics  of  partially  coherent  light  [31] 


4  =  A-(l  +  ^).  (4.40) 

where  K  is  the  distribution  mean  and  M.  is  the  estimated  speckle  parameter.  A  value 
of  Af  =  60  was  estimated  from  the  experimental  data  and  used  for  the  generation  of 
simulated  data  [54],  The  results  of  the  comparison  are  tabulated  in  Table  4.2. 


Table  4.2:  This  table  describes  the  simulated  and  estimated  values  for  Fried’s  seeing 
parameter  tq  as  estimated  using  a  blind  deconvolution  algorithm  that  uses  only  the 
short-exposure  OTF  of  Eqn.  4.1  compared  to  the  same  algorithm  using  the  total 
system  OTF  described  by  Eqn.  4.17. 


Simulated  rg  in  cm 

estimated  rg  using 

estimated  vq  using  AOTF 

5 

4.7 

5.1 

8 

7.5 

7.8 

10 

9.6 

9.8 

15 

14.4 

14.7 

20 

19.1 

19.7 

The  additional  anisoplanatic  blur  components  modeled  by  the  AOTF  increased 
the  accuracy  of  the  estimation  of  Fried’s  parameter  from  5%  mean  error  to  within  2% 
using  simulated  imagery. 

4.2.4  Results  Obtained  using  Experimentally  Colleeted  Image  Data.  Exper¬ 
imentally  collected  data  was  limited  to  hve  pairs  of  300-image  datasets  collected  for  a 
particular  atmospheric  condition  on  a  controlled  mountain-top  optical  range.  A  reso¬ 
lution  bar  target  and  a  step-intensity  target  were  imaged  according  to  the  parameters 
outlined  in  Table  3.1.  The  gated  laser  imaging  camera  was  located  atop  the  North 
Oscura  Peak  site  at  the  White  Sands  Missile  Range,  New  Mexico.  The  site  elevation 
was  7993  feet  (2436  m)  MSL,  while  the  target  site  (Beck  Site)  was  located  at  a  height 
of  5060  ft  (1542  m)  MSL.  The  slant-path  range  to  target  was  10,040  meters.  This 
geometry  resulted  in  a  downlook  angle  of  approximately  17  degrees.  Weather  condi- 
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tions  were  extremely  dry,  with  a  humidity  of  16%.  Individual  images  were  collected 
with  a  gate  time  of  12  ns,  at  a  frame  collection  rate  of  10  Hertz.  Both  targets  were 
arranged  such  that  only  slight  azimuth  change  was  required  to  image  either  target, 
ensuring  similar  atmospheric  prohles.  Additionally,  imaging  of  each  pair  of  target  sets 
was  separated  in  time  by  less  than  three  minutes.  The  10  km  optical  path  to  the  re¬ 
mote  target  prevented  accurate  atmospheric  truth  using  scintillometer  measurements. 
In  order  to  compare  results,  the  atmospheric  seeing  condition  of  the  experimentally 
imaged  step-intensity  target  was  estimated  using  the  knife-edge  OTF  estimation  tech¬ 
nique  described  in  Sec.  4.2. 1.2.  Wind  speed  was  recorded  in  excess  of  35  meters  per 
second  at  the  optical  aperture,  validating  the  assumption  of  independent  turbulence 
realizations  for  each  of  the  10  Hertz  frame  rate  speckle  images. 

Blind  deconvolution  of  motion-compensated  frame  averages  of  the  hve  resolu¬ 
tion  bar  target  data  sets  was  performed  using  the  MAP  algorithm  briefly  described 
in  Sec.  4.2.2.  The  estimated  seeing  conditions  were  compared  for  both  the  short- 
exposure  OTF  and  the  new  AOTF  system  models  for  each  data  set.  These  hgures 
were  compared  to  atmospheric  truth  estimates  derived  from  the  knife-edge  estimation 
technique  outlined  in  Sec.  4.2. 1.2,  and  tabulated  in  Table  4.3. 


Table  4.3:  Table  describing  the  estimated  values  for  Fried’s  seeing  parameter  tq 

for  experimentally  collected  imagery  as  estimated  using  a  knife-edge  OTF  technique, 
as  well  as  a  blind  deconvolution  algorithm  that  uses  only  the  short-exposure  OTF  of 
Eqn.  4.1  compared  to  the  same  algorithm  using  the  total  system  OTF  described  by 
Eqn.  4.17. 


Knife-Edge  tq  in  cm 

estimated  ro  using  Hse 

estimated  ro  using  AOTE 

3.9 

3.6 

3.8 

4.1 

3.4 

3.7 

4.6 

4.3 

4.7 

3.2 

3.0 

3.2 

3.6 

3.4 

3.6 

Fig.  4.5  shows  a  representative  10  km  motion-compensated  frame-average  image 
of  the  resolution  bar  target  as  input  to  the  algorithm,  while  Fig.  4.6  (a)  shows  the 
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Figure  4.5:  Motion-compensated  frame  average  (MCFA)  image  created  from  50 

frames  of  experimentally  collected  laser  radar  data  for  input  to  the  MAP  blind  de- 
convolution  algorithm.  The  seeing  condition  estimated  using  knife-edge  techniques 
was  found  to  be  tq  =  3.9. 

result  of  deconvolution  using  the  short  exposure  transfer  function.  Figure  4.6  (b) 
illustrates  the  same  image  deconvolved  using  the  improved  AOTF  system  model. 

The  additional  anisoplanatic  blur  components  modeled  by  the  AOTF  increased 
the  accuracy  of  the  estimation  of  Fried’s  parameter  from  8.6%  mean  error  to  within 
2.9%  using  experimentally  collected  image  data. 

4-3  Conclusions  and  Discussion 

The  results  presented  in  Tables  4.2  and  4.3  indicate  that  blind  deconvolution 
estimation  of  tq  was  slightly  pessimistic  (poorer  seeing  condition)  when  only  the  short- 
exposure  average  OTF  was  assumed  for  a  valid  imaging  model.  The  inclusion  of  the 
AOTF  into  the  system  OTF  tended  to  yield  less  pessimistic  estimates  of  the  seeing 
parameter  that  were  closer  to  truth  conditions  for  both  the  synthetic  and  experimen¬ 
tally  collected  imagery.  The  addition  of  an  anisoplanatic  related  uncorrelated  blur 
component  to  the  deconvolution  kernel  appears  to  increase  the  accuracy  of  estimated 
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(a)  Isoplanatic  MAP  estimate  image  of  experi-  (b)  Anisoplanatic  MAP  estimate  image  of  exper- 
mentally  collected  laser  radar  image  data  imentally  collected  laser  radar  image  data 


Figure  4.6;  Comparison  of  the  MAP  estimated  image  using  the  isoplanatic  and 
anisoplanatic  deconvolution  kernel  functions.  Subhgure  (a)  is  MAP  estimate  image  of 
experimentally  collected  laser  radar  image  data  using  short-exposure  OTF  only.  Esti¬ 
mated  fo  =  3.6.  Subhgure  (b)  is  the  MAP  estimate  image  of  experimentally  collected 
laser  radar  image  data  using  combined  anisoplanatic  system  model  of  Eqn.  4.17.  Es¬ 
timated  tq  =  3.8.  Note  slightly  increased  image  detail  that  is  apparent  in  the  smaller 
bar  patterns  compared  to  Subhgure  (a). 
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seeing  condition  over  a  fairly  broad  range  of  simulated  seeing  conditions,  as  well  as  at 
least  a  limited  range  of  seeing  conditions  in  the  case  of  experimentally  collected  data. 

The  additional  accuracy  gained  by  improved  estimation  of  atmospheric  condi¬ 
tions  might  be  beneficial  in  many  types  of  imaging  applications.  For  example,  a  novel 
seeing  monitor  might  be  constructed  that  requires  only  the  collection  of  random  im¬ 
ages  with  sufficiently  high  spatial  frequency  content,  rather  than  specihc  test  patterns 
designed  to  assist  the  seeing  monitor.  It  should  be  noted  that  the  derived  AOTF  is  not 
limited  to  coherent  imagery,  and  the  estimation  algorithm  outlined  in  [54]  is  equally 
applicable  to  incoherent  imagery  if  slightly  modihed  to  estimate  Poisson  statistics 
rather  than  the  negative  binomial  statistics  of  partially  coherent  laser  illumination. 
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V.  Weighted  Averaging  of  an  Ensemble  of  Collected  Image 

Frames 


The  research  described  in  previous  chapters  has  introduced  novel  techniques  for 
remote  scene  and  atmospheric  seeing  condition  estimation.  To  use  the  methods 
described,  several  short  duration  images  collected  by  a  suitable  coherent  vision  sys¬ 
tem  must  be  combined  in  order  to  reduce  the  deleterious  effects  of  atmospheric  and 
coherent  speckle.  The  need  to  combine  multiple  frames  of  image  data  introduces  the 
question  of  how  best  to  effect  this  combination.  In  one  extreme,  the  image  processor 
might  choose  to  simply  average  some  number  of  available  images  without  regard  to 
translational  image  alignment  or  registration.  Without  such  registration,  the  compos¬ 
ite  image  may  be  considered  the  resultant  output  of  an  optical  system  well  described 
by  the  long-exposure  OTF  [31]. 

Image  registration  dramatically  increases  high  spatial  frequency  image  detail  in 
the  resulting  average  image  by  reducing  globally  distributed  motion  blur  due  to  cam¬ 
era  platform  vibration  or  tilt  and  tip  components  of  a  turbulent  atmospheric  viewing 
path.  In  the  case  where  viewing  conditions  are  isoplanatic,  the  optical  system  can  be 
described  as  shift  invariant,  and  almost  all  such  motion  blur  can  be  removed  from  the 
averaged  image  given  a  sufficiently  robust  image  registration  algorithm.  The  resulting 
average  image  may  be  considered  the  product  of  a  shift-invariant  optical  system  with 
an  OTF  well  described  using  the  short-exposure  OTF  [31].  However,  there  exist  fun¬ 
damental  as  well  as  practical  limitations  to  the  overall  image  improvement  realized 
by  image  registration  [58] .  Given  such  limitations,  it  is  clear  that  some  level  of  global 
motion  blur  will  remain  in  the  composite  image. 

In  addition  to  the  image  degradation  caused  by  unresolved  global  spatial  regis¬ 
tration,  there  exists  the  problem  of  local  subimage  blur  due  to  the  effects  of  anisopla- 
natism.  Such  effects  were  the  subject  of  study  in  Chapter  IV.  Due  to  the  shift-variant 
system  model  that  describes  individual  image  frame  formation  through  a  wide  FOV 
optical  system  operating  in  a  turbulent  atmospheric  environment,  global  image  regis¬ 
tration  techniques  are  inadequate  to  deal  with  the  isolated  and  relatively  uncorrelated 
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motion  of  the  myriad  of  isoplanatic  patches  that  comprise  each  image  frame.  In  cases 
of  heavy  turbulence,  individual  frames  may  become  severely  warped  or  distorted  due 
to  the  shift-variant  imaging  system.  The  inclusion  of  such  frames  into  the  ensemble 
average  image  will  manifest  as  image  blur,  since  individual  frames  are  statistically 
uncorrelated,  and  the  turbulent  motion  of  the  atmosphere  yields  many  such  frames 
with  distinctly  warped  image  areas. 

This  chapter  supplements  the  research  presented  in  prior  sections  by  seeking 
a  useful  method  whereby  the  effects  of  such  outlier  frames  are  minimized.  A  brief 
background  of  frame  selection  for  optical  systems  and  outlier  detection  in  registration 
algorithms  is  presented  in  Sec.  5.2.  A  cost  function  is  developed  using  maximum 
likelihood  estimation  theory  in  Sec.  5.3  that  quantihes  the  admissibility  of  a  particular 
frame  to  the  overall  ensemble  frame  average.  The  likelihood  function  is  maximized 
using  an  iterative  algorithm  and  compared  to  a  simpler  model  of  the  imaging  process 
that  admits  a  direct  solution  of  the  maximization  problem.  The  research  leads  to 
an  elegant  binary  hypothesis  that  may  be  used  to  discard  frames  from  the  ensemble 
depending  on  a  simple  likelihood  ratio  test.  The  outlier  estimator  performance  is 
demonstrated  using  both  synthetically  generated  as  well  as  experimentally  collected 
image  data  in  Sections  5.5  and  5.6.  The  resulting  development  allows  for  signihcant 
image  enhancement  in  cases  where  atmospheric  distortion  and  image  registration 
dehciencies  cause  degradation  in  the  frame  average  image,  as  discussed  in  Sec.  5.7. 

5. 1  Image  Improvement  by  Averaging  an  Ensemble  of  Registered  Speckle 
Image  Frames 

Considering  the  problem  of  image  reconstruction  from  a  sequential  series  of  short 
time-gated  exposures  from  a  coherent  laser  imaging  system,  accurate  translational 
registration  is  a  critical  initial  requirement  due  to  the  atmospheric  tip  and  tilt  that 
is  characteristic  of  long  viewing  paths  to  the  target  scene,  as  well  as  motion  and 
vibration  encountered  by  the  imaging  platform.  Accurate  frame  registration  permits 
the  subsequent  reduction  of  image  speckle,  a  signihcant  concern  for  coherent  imaging 
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systems  that  use  illuminating  sources  with  long  coherence  times.  In  addition  to  laser 
speckle,  the  atmospheric  turbulence  between  the  target  and  imaging  system  causes 
individual  ensemble  frames  to  become  heavily  corrupted  by  speckle  due  to  the  various 
phase  delays  imparted  by  the  non-uniform  atmospheric  index  of  refraction.  To  further 
complicate  the  problem,  long  distance  viewing  often  results  in  very  low  photon  counts 
despite  employment  of  high-powered  illuminating  laser  sources.  The  summation  of 
several,  well-registered  image  frames  is  often  mandatory  to  increase  image  SNR  and 
reduce  image  speckle. 

The  individual  images  returned  from  a  coherent  imaging  system  over  long  prop¬ 
agation  paths  with  long  illuminator  coherence  times  are  typihed  by  high  levels  of 
speckle  and  large  intensity  variance.  Figure  5.1  shows  a  cropped  128x128  pixel  por¬ 
tion  of  a  single  image  of  a  resolution  bar  target  board  collected  by  an  experimental 
laser-vision  system  at  a  range  of  approximately  10  kilometers.  Typical  imaging  sys¬ 
tems  permit  frame  rates  in  the  tens  of  Hertz,  producing  image  ensembles  ranging 
from  10  to  perhaps  100  images  over  some  acceptably  short  dwell  period.  As  is  evident 
from  Fig  5.1,  speckle  and  intensity  variance  makes  automated  image  registration  dif- 
hcult.  However,  the  method  of  fast  vector  projection  correlation  has  been  applied  to 
this  problem  with  remarkable  success  [7,53].  Figure  5.2  illustrates  the  improvement 
gained  by  automatically  registering  and  averaging  an  ensemble  of  50  image  frames 
collected  under  the  same  conditions  as  shown  in  Fig  5.1. 

Correlation-based  registration  processing  is  not  without  occasional  error,  since 
estimation  of  the  shift  parameters  for  some  frames  in  the  ensemble  may  be  hampered 
by  false  correlation  peaks.  Such  is  often  the  case  when  specular  glint  in  a  partic¬ 
ular  frame  erroneously  correlates  with  actual  bright  features  in  the  remote  scene 
image.  This  effect  may  be  noted  by  the  high-intensity  return  in  the  lower  corner  of 
Fig.  5.2,  where  a  false  correlation  peak  occurred  due  to  a  bright  specular  return  from 
an  off-screen  portion  of  the  image.  In  such  cases,  it  is  important  to  recognize  that  a 
registration  error  has  indeed  occurred,  and  that  steps  are  taken  to  either  re-register 
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Figure  5.1;  Intensity  image  formed  from  a  single  frame  of  reflected  coherent  light 
from  a  resolution  board  target  collected  at  an  experimental  optical  range  from  a 
distance  of  10  kilometers. 

the  frame  or  to  deemphasize  the  contribution  of  the  frame  to  the  ensemble  average 
image. 

5.2  Background 

This  research  uses  a  maximum  likely  (ML)  formulation  to  establish  weighting 
coefficients  for  individual  frames  in  a  multi-frame  image  ensemble.  Several  research 
teams  have  considered  ML  techniques  for  the  development  of  novel  and  robust  image 
registration  algorithms  [11,17,32,39,81,82].  Although  considerable  research  has  been 
conducted  to  develop  new  and  enhanced  image  registration  algorithms,  the  literature 
is  sparse  with  general  techniques  that  quantify  the  goodness-of-£t  of  a  particular  image 
frame  relative  to  the  average  image.  Several  Bayesian  treatments  of  outlier  detection 
within  the  general  context  of  image  processing  are  available  [27,34],  however,  these 
methods  assume  statistical  models  which  are  not  necessarily  applicable  to  the  random 
processes  that  govern  the  detection  of  coherent  illumination  through  a  non-uniform 
atmosphere. 
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Figure  5.2:  Frame  average  of  50  consecutive  image  frames  collected  at  an  experimen¬ 
tal  optical  range  from  a  distance  of  10  kilometers.  Although  substantial  improvement 
is  evident  when  compared  to  Fig.  5.1,  registration  error  may  be  noted  by  the  bright 
specular  return  error  that  is  seen  in  the  lower  left  corner  of  the  target  board.  These 
errors  are  caused  by  false  correlation  peaks  dne  to  a  bright  specular  return  from  the 
door  handle  of  the  supporting  truck  (off-screen  to  left).  Some  vertical  ghosting  of  the 
upper  bars  is  also  evident. 

Fried  realized  that  the  probability  of  obtaining  accnrate,  high  spatial-frequency 
images  from  an  optical  system  decreases  exponentially  as  the  ratio  of  the  optical  aper¬ 
ture  diameter  to  the  seeing  parameter  {D /tq)  increases  [25].  Several  researchers  have 
snbsequently  attempted  to  establish  frame  selection  performance  limits  for  enhance¬ 
ment  of  reconstrncted  imagery  [22,59],  and  they  offer  a  series  of  image  qnality  metrics 
that  may  be  used  to  compare  images  retrieved  from  an  optical  system.  However,  no 
signihcant  research  can  be  found  that  describes  likelihood-based  methods  to  identify 
which  frames  should  be  considered  candidates  for  removal  from  the  ensemble  before 
averaging. 

It  was  shown  in  Sec.  2.6.1  that  the  pixel  intensity  distribntion  along  a  series  of 
temporally  separated  image  frames  gathered  from  a  partially  coherent  illumination 
system  follows  a  Negative  Binomial  (NB)  random  process  [31].  This  distribution  may 
be  understood  to  be  a  more  general  treatment  of  the  Poisson  intensity  process  com- 
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monly  applied  to  the  case  of  incoherent  illumination.  As  demonstrated  in  Sec.  2.6.1, 
the  summation  of  many  frames  of  NB  distributed  intensity  data  result  in  images  that 
are  also  governed  by  a  NB  process  with  increased  speckle  parameter  A4.  Using  this 
physically  based  statistical  model  of  partially  coherent  illumination,  one  may  con¬ 
struct  a  Bayesian  estimator  that  yields  a  likelihood  function  for  the  weight  of  each 
image  in  a  series  of  independently  collected  image  frames  within  a  temporally  contigu¬ 
ous  ensemble.  It  should  be  pointed  out  that  the  coherence  times  of  candidate  laser 
illumination  systems,  while  long  compared  to  incoherent  sources,  is  actually  rather 
short  compared  to  the  frame  rate  of  typical  laser-vision  systems.  Typical  coherence 
times  of  laser  systems  are  measured  in  the  nanosecond  to  microsecond  range,  while 
frame  collection  frame  rates  are  on  the  order  of  10-30  Hertz.  Very  small  changes  in  an¬ 
gle  between  the  target  and  camera  due  to  platform  motion  or  atmospheric  turbulence 
cause  almost  complete  decorrelation  between  images  collected  during  each  frame  pe¬ 
riod.  To  a  very  good  approximation,  individual  frames  gathered  by  such  systems  may 
be  treated  as  statistically  independent  realizations  of  the  underlying  noise  process. 

The  development  of  an  estimator  for  the  relative  weights  that  might  be  assigned 
to  individual  short-exposure  images  is  based  on  the  inherent  distinction  between  image 
degradation  that  occurs  due  to  the  coherent  imaging  process  versus  that  caused  by 
anisoplanatic  image  warping  or  poor  image  frame  registration.  Although  difficult  to 
identify  due  to  the  lack  of  a  priori  knowledge  of  the  remote  scene  and  atmospheric 
seeing  conditions,  it  may  be  assumed  that  the  intensity  distribution  of  individual 
pixels  due  to  image  warping  or  mis-registration  does  not  follow  a  NB  distribution.  A 
possible  exception  to  this  assumption  might  occur  for  the  case  of  imaging  a  scene  of 
uniform  reflectance.  However,  for  the  vast  majority  of  interesting  cases,  there  exists 
no  mechanism  to  cause  one  to  believe  that  such  a  distribution  would  be  governed 
by  a  NB  noise  process  in  the  more  general  case.  Under  this  framework,  it  is  easy 
to  understand  why  the  likelihood  of  an  individual  image  might  be  assigned  a  lower 
value  in  those  cases  where  warping  or  mis-registration  has  occurred.  If  no  such  image 
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degradation  exists,  a  candidate  frame  will  be  well  modeled  by  the  statistics  of  the  NB 
process,  and  the  assigned  likelihood  will  be  relatively  high. 


5.3  Frame  Weight  Estimator  Development 

A  series  of  J  images  is  collected  by  a  system  that  propagates  highly  coherent 
light  towards  a  target  scene  through  a  volume  of  turbulent  atmosphere  and  reflects 
from  the  target  back  to  the  optical  receiver  aperture  as  discussed  in  Sec.  2.1.  A 
physically  motivated  model  for  the  statistical  variation  of  pixel  intensity  measured  in 
photons  is  negative  binomial  [31], 
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where  K  is  the  mean  photon  count,  M.  is  the  speckle  parameter,  and  K  is  the  random 
photon  count  at  the  detector.  Let  K  =  dj{x  —  aj,y  —  (3j),  the  number  of  photons 
arriving  at  the  detector  for  each  pixel  for  the  detected  image  in  the  ensemble. 
Here,  oij  and  jdj  are  the  previously  estimated  shifts  for  each  of  the  J  image  frames  in 
the  ensemble  according  to  some  arbitrary  registration  algorithm. 


In  order  to  incorporate  the  effects  of  frame  weighting  within  an  ensemble  of 
images,  the  mean  pixel  intensity  may  be  modeled  as  the  weighted  average  of  each  of 
the  ensemble  images  shifted  according  to  the  estimated  registration  components,  cxj 
and  jSj.  Let  K  =  i  {x,y),  dehned  mathematically  as 
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where  A  is  the  vector  of  J  weights  that  remain  to  be  estimated  according  to  the 
statistics  of  Eqn.  5.1,  and  n  is  used  to  index  the  sum  to  avoid  confusion  later  in 
this  derivation.  Thus,  i{x,y)  may  be  thought  of  as  a  pixel  of  the  weighted  motion- 
compensated  frame  average  image.  According  to  this  model,  frames  that  are  found 
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to  have  low  estimates  for  the  corresponding  weight  Aj  will  tend  to  contribute  less  to 
the  frame- averaged  image. 


5.3.1  Maximum  Likelihood  Frame  Weight  Estimator  Derivation.  Without 
the  beneht  of  a  priori  information  on  the  frame  weights,  a  MAP  estimator  cannot  be 
constructed,  however,  a  suitable  maximum  likelihood  (ML)  estimator  may  be  derived. 
Bayes  rule  may  be  used  to  maximize  the  probability  of  frame  weights  given  the  image 
data,  Pr  [Ajidj  {x^y)] ,  by  simply  maximizing  the  probability  of  the  image  data  given 
the  weights,  Pr  [dj  (x,  y)  \Aj\  [71].  The  pixels  of  the  data  are  assumed  independent  and 
identically  distributed  as  is  often  the  case  with  similar  developments  [39],  allowing 
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With  the  assumption  that  the  frames  are  independent  due  to  the  relatively  long 
interframe  period  and  therefore  contain  uncorrelated  speckle,  the  total  probability 
may  be  found  by  multiplying  the  probabilities  of  all  J  image  frames  in  the  ensemble. 
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It  is  convenient  to  maximize  the  natural  logarithm  of  Eqn.  5.4  due  to  numerical 
difficulties  encountered  while  multiplying  many  large  numbers.  The  resulting  log- 
likelihood  function  L  (d)  =  InP^iA  (D|A)  is 
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The  intended  maximization  of  the  log-likelihood  is  with  respect  to  the  weights 
of  the  individual  frames  in  the  ensemble.  Let  Aj^  be  an  arbitrary  frame  weight  in  the 
set  of  Aj.  Note  that  the  first  term  within  the  braces  of  Eqn.  5.5  bears  no  dependence 
on  the  frame  weights  and  may  be  disregarded  in  the  maximization  analysis.  However, 
the  remaining  two  terms  have  an  implicit  relationship  with  Aj  as  dehned  by  Eqn.  5.2. 
After  some  arithmetic,  the  derivative  of  the  log-likelihood  with  respect  to  an  arbitrary 
weight  in  the  ensemble  can  be  expressed  as 
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In  order  to  simplify  this  expression,  it  is  necessary  to  examine  the  change  in 
i{x,y)  with  respect  to  Aj^.  By  substitution  with  Eqn.  5.2, 
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the  change  in  i{x,y)  with  respect  to  Ajg  reduces  to  the  single  term  after  expanding 
the  sum, 

y)  =  {x,  y)  =  Ijdj^  {x  -  aj„y  -  pjP  .  (5.8) 

After  a  bit  more  arithmetic,  the  derivative  may  be  expressed  as 
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By  rearranging  the  terms  and  setting  to  zero  to  find  the  maximum  of  the  log- 
likelihood,  one  obtains 
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Although  an  explicit  solution  for  the  elements  of  A  does  not  appear  to  be  readily 
obtainable,  an  iterative  solution  may  be  found  as  follows.  For  each  particular  frame 
weight  in  the  ensemble,  Ajo,  a  new  estimate,  A”®"'  may  be  recursively  found  from  the 
previous  estimate,  Ajj^'^,  given  a  particular  dataframe  from  the  ensemble,  djo- 

The  update  equation  may  be  conveniently  expressed  as 
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where 

1 

(x,  y)  =  J  X]  A°J^‘^dn  -an,y-  (3^  ■ 

n=\ 

This  recursive  technique  has  the  benefit  of  restricting  A”®"'  to  be  a  member  of 
the  set  of  positive  real  numbers.  To  estimate  the  J  frame  weights  of  images  in  an 
ensemble,  Eqn.  5.11  must  be  iterated  for  each  frame  under  consideration.  At  each 
iteration,  the  average  image  i{x,  y)  is  reformed  with  the  frame  weights  found  from  the 
previous  iteration  according  to  Eqn.  5.2. 
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Appendix  A  describes  a  direct  solution  of  a  likelihood  function  based  on  simpler 
Gaussian  statistics  rather  than  negative  binomial  statistics. 

5.3.2  ML  Frame  Weight  Estimator  Implementation.  The  iterative  algorithm 
was  implemented  using  a  general-purpose  computer.  Note  that  the  numerator  and 
denominator  of  the  update  equation  in  Eqn.  5.11  result  in  strictly  positive  scalar 
quantities  due  to  the  summation  over  ah  frames  and  pixels  in  the  resulting  images. 
A  nominal  starting  point  for  the  frame  weights  is  a  vector  of  ones.  In  practice,  the 
resulting  frame  weights  may  be  normalized  to  sum  to  J  after  each  iteration,  such  that 
likely  frames  within  the  ensemble  remain  close  to  unity,  while  frames  that  are  deemed 
outliers  tend  to  drop  to  values  less  than  unity  as  the  iterations  progress. 

For  all  simulated  and  experimentally  collected  datasets  analyzed,  the  algorithm 
appeared  to  slowly  arrive  at  a  solution  where  outlier  frames  have  very  small  associated 
weights,  while  frames  that  ht  well  to  the  ensemble  mean  have  associated  weights  that 
remain  close  to  unity.  However,  no  clear  strategy  for  terminating  the  iterative  process 
was  discovered.  Prior  to  using  the  algorithm,  an  estimate  for  the  speckle  parameter 
of  the  system,  Af,  must  be  found.  If  enough  frames  of  registered  image  data  are 
available  from  a  calibration  dataset,  M.  may  be  found  directly  from  the  variance 
of  pixel  intensity  along  columns  of  bright  pixels  according  to  the  expression  for  the 
variance  of  the  negative  binomial  distribution  [31], 

Alternatively,  the  maximum  likelihood  estimation  procedure  for  the  speckle  parameter 
described  in  Sec.  2.5  may  be  used. 

Algorithm  convergence  was  notably  faster  for  larger  values  of  the  coherent  source 
speckle  parameter.  Although  the  convergence  of  the  iterative  algorithm  can  be  slow 
for  low  values  of  AI,  an  important  observation  was  that  the  gradients  of  the  weights 
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could  be  inspected  after  only  a  few  iterations  in  order  to  decide  if  the  associated 
frames  indicated  poor  ensemble  registration. 

As  an  example,  Fig.  5.3  shows  the  hrst  difference  of  the  weight  vector  after  only 
two  iterations  using  imagery  collected  from  an  experimental  optical  range.  These  are 
the  same  data  used  to  create  the  images  shown  in  Figures  5.1  and  5.2.  Note  that  hve 
of  the  frames  are  associated  with  weights  that  have  highly  negative  gradients.  These 
frames  correspond  to  the  registration  errors  that  occur  due  to  false  correlation  of  a 
bright  specular  return  source  on  the  remote  scene.  As  the  algorithm  is  allowed  to 
progress  until  the  likelihood  equation  approaches  an  asymptotic  value,  the  weights  of 
the  corresponding  negative  gradients  eventually  approach  zero.  However,  the  num¬ 
ber  of  iterations  required  to  reach  such  convergence  might  be  prohibitive  in  many 
applications.  In  such  applications,  a  simple  method  might  be  devised  to  compare  the 
gradient  of  the  weights  to  the  ensemble  mean  or  median  of  the  gradients  after  several 
iterations,  effectively  detecting  frames  that  do  not  £t  well  to  the  image  ensemble. 
Figure  5.4  (a)  shows  a  typical  weight  corresponding  to  an  outlier  frame  as  a  function 
of  the  number  of  iterations,  while  Fig.  5.4  (b)  shows  the  ensemble  likelihood  as  calcu¬ 
lated  for  each  iteration  step  using  Eqn.  5.5.  Figure  5.5  (b)  shows  the  resulting  frame 
average  image  with  weights  derived  from  the  iterative  algorithm  after  1200  iterations. 
Notably  absent  is  the  contribution  from  the  poorly  registered  frames  visible  in  Fig.  5.2 
and  repeated  for  comparison  in  Fig.  5.5  (a). 

The  outlier  detection  algorithm  did  not  appear  to  be  overly  sensitive  to  poor 
estimates  of  the  speckle  parameter.  While  low  estimates  of  the  system  speckle  param¬ 
eter  did  tend  to  slow  the  rate  at  which  frame  weights  decreased  below  unity,  higher 
estimates  tended  to  speed  weight  decrease  by  a  commensurate  amount.  However,  in 
all  cases  studied,  the  relative  gradient  of  the  weights  seemed  to  provide  a  reliable  and 
robust  indicator  of  frames  that  were  either  poorly  registered  or  had  significant  image 
warping. 
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Figure  5.3:  Plot  of  the  gradient  of  the  weights  (first  difference)  calculated  after  two 
iterations  of  Eqn.  5.11  for  each  of  the  50  weights.  Frames  24,  26,  28,  38  and  48  are 
good  candidates  for  elimination,  dne  to  a  snbstantial  negative  trend  of  the  weights 
towards  zero.  Image  data  was  experimentally  collected  at  a  range  of  10  km  and  had 
an  estimated  speckle  parameter  oi  M.  =  60.  The  vector-correlation  method  of  [7]  was 
nsed  to  register  the  imagery. 

5.4  Frame  Average  Image  Improvement  by  Discarding  Suspect  Outlier 
Image  Frames 

The  frame  weight  estimation  techniqne  described  in  Sec.  5.3  may  be  impracti¬ 
cal  for  some  applications  dne  to  several  shortcomings.  The  most  notable  detractor 
of  the  algorithm  is  its  compntational  burden.  Although  the  component  mathemati¬ 
cal  operations  are  simple,  the  nnmber  of  operations  increases  dramatically  for  large, 
operationally  representative  imagery,  and  large  nnmbers  of  images  within  each  en¬ 
semble.  Another  detractor  involves  the  means  to  stop  snch  an  iterative  algorithm. 
Although  one  might  find  a  suitable  stopping  criteria  by  inspection  of  the  iterative 
weight  differences  or  perhaps  the  rate  of  ascent  of  the  log-likelihood  curve,  no  mea- 
snre  of  optimality  is  guaranteed  by  such  ad-hoc  criteria.  To  overcome  such  limitations, 
the  problem  was  recast  as  a  binary  hypothesis  as  discussed  below. 
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Ensemble  Log-Likelihood  as  a  tunction  of  Iteration  Number 


(a)  Typical  outlier  continuous  frame  weight  as  a 
function  of  iteration. 


(b)  Overall  ensemble  likelihood  as  a  function  of 
iteration. 


Figure  5.4:  Calculated  weight  of  a  typical  outlier  frame  and  ensemble  likelihood  as 
a  function  of  the  number  of  iterations.  As  iterations  progress,  the  weight  assigned 
to  an  outlier  frame  tends  to  decrease  asymptotically  to  zero  (a),  while  the  overall 
ensemble  likelihood  tends  to  increase  (b). 


Signihcant  insight  may  be  obtained  by  examination  of  frame  weight  values  after 
a  large  number  of  estimator  iterations.  For  datasets  which  contain  outliers  due  either 
to  local  or  global  mis-registration,  the  weights  corresponding  to  outlier  frames  tended 
to  decrease  to  very  small  values.  This  observation  suggests  the  utility  of  an  algorithm 
that  permits  only  binary  frame  weighting.  Such  a  system  would  assign  a  frame  weight 
of  unity  for  those  frames  that  £t  well  to  the  ensemble  average,  while  assigning  a  value 
of  zero  to  those  that  did  not. 

The  problem  may  be  simplified  by  formulating  two  distinct  hypotheses.  Let 
represent  the  hypothesis  that  the  frame  is  well  matched  to  the  ensemble  average, 
while  represents  the  case  that  the  frame  is  an  ontlier  relative  to  other  images 
within  the  ensemble.  Under  this  framework,  two  distinct  frame  average  images  may 
be  constructed.  Under  the  candidate  frame  Aj^  should  be  included  in  the  ensem¬ 
ble  and  be  assigned  full  weighting  of  Aj^  =  1.  In  this  case  the  frame  average  may 
be  constructed  by  including  the  frame  in  the  set  of  equally  weighted  frames  in  the 
ensemble. 
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(a)  Unweighted  average  image.  (b)  ML  frame  weighted  average  image. 

Figure  5.5:  Comparison  of  the  weighted  average  image  created  by  applying  the 

iteratively  determined  weights  to  the  frames  within  the  ensemble,  to  the  unweighted 
average  image  using  all  frames  in  the  ensemble.  Subimage  (a)  shows  the  unweighted 
average  image.  Due  to  the  low  value  of  some  of  the  weights,  the  resulting  weighted 
average  image  of  subimage  (b)  is  essentially  created  by  eliminating  those  weights  that 
have  driven  to  values  close  to  zero  after  1200  iterations.  Note  the  absence  of  the 
bright  specular  return  and  vertical  ghosting  of  horizontal  resolution  bars  evident  in 
subimage  (a). 
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For  the  case  of  the  alternate  hypothesis,  the  frame  average  may  be  con¬ 
structed  by  deleting  the  frame  under  consideration,  while  assigning  equal  weighting 
to  the  remaining  frames  in  the  ensemble. 
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To  compute  a  binary  hypothesis  test,  the  probability  distributions  governing 
each  hypothesis  must  be  computed  and  compared.  The  likelihood  ratio  test  provides 
a  convenient  method  to  effect  such  a  comparison  [71].  Under  this  construct,  if 
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where  7  represents  some  threshold  determined  by  cost,  then  the  data  was  most  likely 
generated  under  hypothesis  H^.  The  probability  distributions  Pr[H^]  and  Pr[H^] 
may  be  quantihed  by  substitution  of  Equations  5.13  and  5.14  into  Eqn.  5.1  with  K 
equal  to  either  (x,  y)  or  U  (x,  y).  For  the  frame  of  the  ensemble,  the  correspond¬ 
ing  distributions  are 
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for  the  the  and  hypotheses  respectively. 

Upon  assumption  of  statistically  independent  ensemble  images,  the  total  prob¬ 
ability  likelihood  ratio  for  the  frame  may  be  written  as 
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Although  not  identical,  the  leading  gamma  terms  in  both  the  numerator  and 
denominator  may  be  simplihed  by  noting  that  all  but  the  jg^  term  will  survive  the 
division  operation  of  the  likelihood  ratio.  The  resulting  expression  may  then  be 
written  as 
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is  the  surviving  gamma  term  after  division.  To  avoid  numerical  difficulties,  a  log- 
likelihood  ratio  may  be  calculated  as 
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To  decide  if  the  suspect  frame  should  be  retained  for  inclusion  within  the 
weighted  frame  average,  one  needs  only  compute  the  ratio  for  each  of  the  J 
frames  in  the  ensemble.  As  will  be  clarihed  in  the  following  section,  inspection  of 
the  distribution  of  likelihood  ratios  reveals  those  frames  that  should  be  retained  or 
discarded  in  order  to  increase  the  overall  likelihood  of  the  ensemble  averaged  image. 

5.4-1  Distribution  of  Likelihood  Ratios  for  Ensemble  Images.  The  random 
nature  of  the  images  collected  by  the  coherent  imaging  system  yields  a  likelihood 
ratio  that  may  also  be  considered  a  random  variable  D.  In  order  to  decide  which 
frames  should  be  discarded,  a  rule  must  be  established  to  compare  the  elements  of 
the  likelihood  ratio  ©  to  some  threshold  7.  A  simple  rule  might  set  the  threshold  to 
the  sample  mean  of  the  likelihood  ratio  population,  and  discard  those  frames  that  fall 
below  this  mean, 

(5'23) 

Using  such  a  rule,  if  ©jq  >  D  then  the  hypothesis  must  be  declared  as  indicated 
in  Eqn.  5.15  for  the  frame. 

Unfortunately,  such  a  simple  rule  presents  difficulties  in  practical  applications. 
An  important  shortcoming  of  this  simple  rule  may  be  understood  by  considering 
the  case  where  an  ensemble  contains  many  images  that  £t  fairly  well  to  the  average 
ensemble  image,  several  images  that  are  moderately  degraded  by  warping  or  mis¬ 
registration,  and  yet  a  few  images  that  are  severely  degraded.  Using  the  rule  of 


5-18 


Eqn.  5.23,  the  calculated  threshold  would  be  heavily  influenced  by  the  few  severely 
corrupted  image  frames.  In  practice,  this  threshold  might  be  so  large  that  the  mod¬ 
erately  corrupted  image  frames  would  not  be  identified  as  outliers  relative  to  the 
majority  of  the  ensemble  family. 

The  development  of  a  more  useful  decision  rule  requires  an  understanding  of 
the  distribution  of  the  random  variable  D.  Unfortunately,  complete  characterization 
of  this  distribution  requires  a  priori  information  concerning  the  remote  target  scene 
as  well  as  the  current  atmospheric  conditions  under  which  the  imagery  was  collected. 
However,  the  problem  may  be  simplified  by  considering  two  distinct  cases.  One  case 
considers  the  vector  of  likelihood  ratios  D°  corresponding  to  frames  from  an  ensemble 
of  images  corrupted  only  by  the  negative  binomial  noise  process.  The  second  case 
involves  likelihood  ratios  corresponding  to  frames  from  an  ensemble  where  some 
image  frames  are  corrupted  by  other  noise  processes  such  as  global  mis-registration 
or  anisoplanatic  frame  warping. 

Considering  the  case  of  an  argument  may  be  made  to  demonstrate  that  the 
distribution  approaches  Gaussian  in  the  limit  as  the  number  of  frames  in  the  ensemble 
grows  large.  The  denominator  of  Eqn.  5.22  may  be  considered  to  be  a  constant  K  for 
a  given  ensemble  of  images,  and  does  not  change  as  a  function  of  the  selected  frame 
jo  under  consideration.  In  this  case,  the  random  variable  may  be  expressed  as 
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In  the  limit  as  the  number  of  statistically  independent  frames  grows  large,  the  summa¬ 
tion  of  J  unknown  distributions  approaches  Gaussian  per  the  Central  Limit  Theorem. 
Division  of  this  random  variable  by  the  constant  K  only  scales  this  Gaussian  distri¬ 
bution.  This  result  is  intuitively  satisfying,  as  it  implies  that  a  subject  frame  under 
the  case  of  D°  will  result  in  a  likelihood  ratio  that  symmetrically  falls  on  either  side  of 
the  mean  of  some  unimodal  distribution.  This  distribution  is  depicted  in  Fig.  5.6  as 
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Pci{D^).  In  some  cases  the  frame  may  fit  better  to  the  ensemble  average,  in  other  cases 
worse.  Such  behavior  may  be  attributed  to  the  unbiased  £t  that  a  particular  frame 
will  have  within  an  ensemble  where  the  only  distortion  is  caused  by  the  statistics 
attributed  to  partially  coherent  illumination. 

In  the  case  of  D^,  frames  that  do  not  fit  well  to  the  ensemble  will  cause  the  nu¬ 
merator  of  Eqn.  5.22  to  become  large,  since  the  likelihood  P^iD)  will  increase  for  such 
frames.  However,  those  frames  that  do  fit  the  average  ensemble  image  will  distribute 
as  for  the  case  of  D°.  This  mechanism  destroys  the  symmetry  of  the  distribution. 
Frames  that  do  not  fit  the  ensemble  will  tend  to  skew  the  distribution  by  pushing 
the  mean  of  the  distribution  towards  larger  values.  This  distribution  is  depicted  in 
Fig.  5.6  as  pd{D^).  This  result  is  useful  to  help  determine  how  best  to  select  those 
frames  that  must  be  discarded  from  the  ensemble. 


Likelihood  Ratios  Corresponding  to  Outliers 


Figure  5.6:  Plot  of  two  distinct  distributions  of  the  random  variable  D.  The  sym¬ 
metrical  normal  distribution  pd{D^)  results  from  the  likelihood  ratio  test  for  the  case 
where  all  frames  are  corrupted  only  by  negative  binomial  noise.  The  skewed  dis¬ 
tribution  corresponding  to  outlier  frames  is  labeled  as  pd{D^).  For  the  case  where 
a  significant  number  of  frames  are  corrupted  by  other  noise  processes  such  as  im¬ 
age  warping  or  global  mis-registration,  the  symmetry  of  the  PDF  pd{D^)  is  destroyed 
because  the  value  of  the  likelihood  ratio  of  Fqn.  5.22  tends  to  increase  for  such  frames. 
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In  contrast  to  the  rule  proposed  in  Eqn.  5.23,  a  more  appropriate  algorithm 
evaluates  and  compares  the  distribution  of  an  unknown  likelihood  ratio  vector  © 
to  a  Gaussian  distribution.  If  the  distribution  appears  Gaussian  within  some  pre- 
dehned  conhdence  metric,  the  algorithm  may  be  terminated  and  only  image  data  with 
likelihood  ratios  greater  than  a  preset  threshold  would  be  discarded.  Alternatively,  in 
the  case  where  the  distribution  of  D  is  found  to  be  sufficiently  distinct  (skewed)  from 
a  Gaussian  distribution,  frames  with  likelihood  ratios  above  the  threshold  should  be 
discarded  and  the  likelihood  ratio  test  repeated.  A  suitable  threshold  may  be  found 
by  Erst  calculating  the  unbiased  sample  variance 

1  ^ 

=  ivVT  E  -  ®)  •  (5-25) 

i=i 

where  D  is  defined  by  Eqn.  5.23.  The  standard  deviation  of  the  distribution  is  ap 
and  may  be  used  to  formulate  a  suitable  threshold.  Setting  a  1-sigma  threshold  of 
7  =  ID)  -|-  (Tp  identifies  those  frames  that  are  reasonably  distant  from  the  process  mean 
for  elimination  from  the  ensemble  average.  After  elimination,  another  ratio  test  must 
be  performed  to  ensure  that  elimination  of  these  outliers  yields  a  Gaussian  distribution 
of  likelihood  ratios.  If  not,  the  algorithm  must  be  repeated  until  outliers  have  been 
eliminated.  Such  processing  avoids  the  possibility  of  undetected  moderate  outliers 
due  to  the  presence  of  frames  that  are  far  removed  from  the  ensemble  average. 

5.4.2  Testing  the  Likelihood  Ratio  for  Gaussian  Distribution.  Several  statis¬ 
tical  tools  exist  to  test  an  unknown  distribution  for  Gaussian  £t.  Of  particular  merit 
is  the  Lilliefors  statistical  test  [16],  which  does  not  require  a  priori  knowledge  of  the 
parameters  of  the  Gaussian  distribution  to  which  the  data  are  compared. 

The  Lilliefors  test  requires  calculation  of  the  sample  mean  and  variance  as  de¬ 
scribed  in  Eqn.  5.23  and  5.25.  The  observed  data  ratios  are  statistically  normalized 
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by  subtracting  the  mean  and  dividing  by  the  variance, 


Zj  — 


(^D 


(5.26) 


This  data  normalization  process  distinguishes  the  test  as  a  rehnement  of  the  Kolmogorov- 
Smirnov  test  for  normality  [16].  Two  hypotheses  are  proposed.  Tio  is  the  hypothesis 
that  the  random  sample  comes  from  a  normal  distribution  with  unspecihed  mean  and 
variance,  while  Tii  denotes  the  hypothesis  that  the  data  comes  from  a  non-normal 
distribution.  The  empirical  cumulative  distribution  S{i)  function  is  constructed  from 
the  normalized  data  samples  and  compared  to  the  standard  normal  cumulative  distri¬ 
bution  function  F{j).  The  Lilliefors  test  statistic  T  is  dehned  as  the  largest  difference 
between  the  two  cumulative  distributions. 


T  =  sup|F(j)-^(j)l-  (5.27) 

j 

Using  this  construct,  reject  Tio  (delcare  the  sample  data  as  non-normally  dis¬ 
tributed)  with  signihcance  level  a  if  the  test  statistic  is  greater  than  p  =  1  —  a.  Tables 
of  Lilliefors  quantiles  for  varying  sample  sizes  and  a  may  be  found  in  the  literature, 
e.g.  [16].  For  sample  sizes  of  J  >  30  samples,  a  p- value  of  95%  can  be  found  by  cal¬ 
culating  tc.gs  =  0.866/y J.  If  T  is  found  to  exceed  tc.gs,  the  hypothesis  that  the  data 
were  drawn  from  a  normal  distribution  may  be  rejected  within  a  conhdence  interval 
of  95%. 


5.5  Results  using  Simulated  Anisoplanatie  Imagery  Data 

Simulation  of  a  coherent  laser  vision  system  was  discussed  in  detail  in  Sec¬ 
tions  2.1  through  2.3.  The  anisoplanatie  viewing  conditions  encountered  by  a  wide 
FOV  optical  system  may  cause  significant  anisoplanatie  image  warping  due  to  the 
shift-variance  of  the  optical  system.  In  addition,  the  speckle  caused  by  coherent  il¬ 
lumination  and  atmospheric  phase  delay  makes  accurate  global  image  registration 
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difficult,  especially  where  long  optical  paths  reduce  the  available  signal-to- noise  ratio 
at  the  imaging  detector.  The  effects  may  be  accurately  simulated  using  the  same  tech¬ 
niques  used  to  conduct  the  remote  scene  and  seeing  condition  estimation  experiments 
documented  in  Chapters  III  and  IV. 

A  resolution  board  target  intensity  pattern  was  propagated  to  the  optical  aper¬ 
ture  through  various  levels  of  turbulence  simulated  by  the  placement  of  Gaussian 
random  thin  phase  screens  along  the  optical  path.  The  synthetically  generated  im¬ 
agery  was  generated  with  a  speckle  parameter  of  Ad  =  60  to  match  the  experimentally 
collected  image  data  discussed  in  Sec.  5.6.  Observation  of  the  individually  generated 
images  under  poor  seeing  conditions  revealed  significant  anisoplanatic  image  warp¬ 
ing  that  was  impossible  to  remove  by  global  image  registration  techniques,  as  may 
be  illustrated  by  the  example  of  Fig.  5.7.  The  vector  correlation  image  registration 
algorithm  described  in  [7]  was  selected  to  remove  global  image  tilt.  Despite  global  tilt 
removal,  signihcant  image  blurring  was  noted  in  the  unweighted  frame  average  due  to 
many  highly  distorted  speckle  images  caused  by  the  random  nature  of  the  generated 
phase  screens.  Figure  5.7  shows  a  comparison  between  a  nominal  simulated  speckle 
image  that  does  not  suffer  dramatic  distortion  effects,  and  an  image  selected  due  to 
its  heavily  distorted  appearance  caused  by  a  particularly  unlucky  phase  screen  that 
was  generated  close  to  the  optical  aperture.  Clearly,  the  removal  of  rogue  images 
such  as  that  shown  in  Fig.  5.7  (b)  from  the  ensemble  average  will  tend  to  enhance  the 
effective  seeing  condition  of  the  system. 

To  understand  the  performance  of  the  detection  algorithm  under  a  broad  range 
of  atmospheric  conditions,  the  simulated  imagery  data  of  Chapter  IV  were  processed 
for  outlier  detection.  The  data  were  partitioned  into  subsets  of  50-frame  image  en¬ 
sembles  for  a  total  of  10  ensembles  per  atmospheric  condition.  D/tq  conditions  of 
10,  4,  5,  2.5,  2,  1.3  and  1  were  simulated  using  a  20  cm  optical  aperture.  Simulation 
parameters  of  these  data  are  summarized  in  Table  5.1.  The  data  were  processed  for 
outliers  using  the  algorithm  detailed  in  Sec.  5.4.  A  threshold  of  one  standard  devia¬ 
tion  above  the  process  mean  was  selected  for  rejection  of  suspected  outliers,  and  the 
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(a)  Nominal  synthetic  image  frame.  (b)  Synthetic  image  heavily  distorted  by  turbu¬ 

lence. 

Figure  5.7:  Comparison  of  a  pair  of  synthetically  generated  resolution  bar  target 

speckle  images  drawn  from  a  randomly  generated  turbulence  simulation.  Simulated 
average  tq  is  2  cm.  over  a  10  kilometer  path.  Subhgure  (a)  shows  a  nominal  frame  from 
the  ensemble.  Distinct  bars  are  apparent  for  all  but  the  smallest  pattern,  however 
some  anisoplanatic  warping  is  notable.  Subhgure  (b)  shows  an  image  frame  heavily 
distorted  by  the  atmospheric  distortion.  Such  a  frame  is  clearly  a  good  candidate  for 
automatic  removal  by  an  outlier  detection  algorithm. 
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process  was  repeated  until  the  likelihood  ratio  appeared  to  be  normally  distributed  as 
discussed  in  Sec.  5.4.1.  In  most  cases  only  two  iterations  of  the  likelihood  ratio  test 
were  required  to  reach  a  normally  distributed  ratio  distribution.  For  good  seeing  con¬ 
ditions,  the  likelihood  ratio  data  appeared  Gaussian  after  the  first  iteration,  and  no 
additional  iterations  were  performed.  In  several  cases  under  poor  seeing  conditions,  a 
third  iteration  was  required  due  to  very  distant  outliers  that  tended  to  mask  outliers 
closer  to  the  process  mean. 


Table  5.1:  Table  describing  the  simulation  parameters  used  to  create  statistically 

accurate  random  speckle  imagery  for  evaluation  of  the  registration  outlier  detection 
algorithm. 


Parameter 

Value 

Slant  Range  to  Target 

10  Kilometers 

Optical  Diameter 

20  Centimeters 

Speckle  Parameter  of  Source 

60 

Pixels  per  Image 

128  by  128 

Pixel  Pitch  of  Detector 

11.8  micrometers 

Mean  Wavelength 

1.54  micrometers 

Focal  Length 

3  Meters 

Images  per  Frame  Ensemble 

50 

Size  of  Imaged  Target  Area 

5  by  5  Meters 

Figure  5.8  depicts  a  typical  improvement  in  the  average  image  gained  by  binary 
frame  weighting  of  the  component  speckle  image  frames.  Intensity  bars  are  included  in 
the  hgures  to  demonstrate  the  additional  contrast  that  is  gained  by  automated  deletion 
of  frames  that  serve  only  to  blur  the  average  image.  Table  5.2  details  the  performance 
of  the  detection  routine  using  simulated  imagery  for  various  atmospheric  conditions. 
The  knife-edge  OTF  estimation  technique  described  in  Sec.  2.9  was  performed  on  the 
weighted  average  images  to  ascertain  the  level  of  improvement  in  spatial  resolution 
gained  by  the  automatic  deletion  of  suspect  frames. 
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Weighted  Average  image 
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(a)  Synthetic  unweighted  average  image.  (b)  Synthetic  binary  weighted  average  image. 

Figure  5.8:  Comparison  of  the  unweighted  (subimage  a)  and  binary  weighted 

(subimage  b)  average  images  for  a  typical  outlier  detection  simulation.  Sample  image 
frames  of  this  50-frame  ensemble  are  shown  in  Figures  5.7(a)  and  5.7(b).  Note  the  in¬ 
creased  image  contrast  indicated  by  the  maximum  photon  count  on  the  intensity  bars, 
as  well  as  the  slight  but  notable  reduction  of  blur  in  the  weighted  average  compared  to 
the  unweighted  average.  Simulated  tq  is  4  cm,  and  21  of  50  frames  were  automatically 
removed  from  the  ensemble  prior  to  computing  the  average  image  shown  in  subimage 
(b). 


Unweighted  Average  image 
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Table  5.2:  Table  showing  the  results  of  weighted  frame  averaging  5  sets  of  syn¬ 

thetically  generated  resolution  bar  imagery  over  a  broad  range  of  atmospheric  seeing 
conditions.  Values  for  number  of  frames  rejected  and  effective  tq  are  averaged  over 
the  ten  50  frame  ensembles  within  each  500  image  data  set. 


Target  Set 

Frames  Rejected 

Actual  To 

Effective  tq 

To  Increase 

Set  1 

24.4 

2  cm. 

2.7  cm. 

35% 

Set  2 

23.4 

4  cm. 

4.7  cm. 

17% 

Set  3 

21.4 

5  cm. 

5.7  cm. 

14% 

Set  4 

18.3 

8  cm. 

8.4  cm. 

5% 

Set  5 

17.3 

10  cm. 

10.2  cm. 

2% 

Set  6 

16.4 

15  cm. 

15.1  cm. 

0.7% 

Set  7 

16.1 

20  cm. 

20.0  cm. 

0% 

5.6  Results  using  Experimentally  Collected  Imagery  Data 

A  large  variety  of  experimentally  collected  image  data  was  available  to  test 
the  outlier  detection  algorithm.  In  addition  to  the  supported  resolution  bar  target 
analyzed  in  Chapters  III  and  IV,  other  interesting  scenes  were  processed  to  test  the 
outlier  detection  algorithm  performance  using  nominal  tactical  imagery.  In  all  cases, 
these  image  datasets  were  pre-processed  for  global  tip  and  tilt  removal  by  the  fast 
vector  correlation  algorithm  described  in  [7].  Figure  5.1  shows  a  nominal  speckle 
image  of  the  supported  resolution  bar  target  observed  at  a  range  of  10  kilometers 
using  a  3  meter  focal  length  optical  system.  Due  to  the  presence  of  large  resolution 
bars,  the  knife-edge  OTF  estimation  technique  describe  in  Sec.  2.9  was  useful  in 
determining  the  amount  of  improvement  that  removal  of  suspect  speckle  image  frames 
produced  for  these  datasets.  Knife-edge  estimation  was  not  available  for  tactical  target 
scenes.  In  these  cases,  the  anisoplanatic  OTF  blind  estimation  algorithms  described 
in  Chapters  III  and  IV  was  used  to  judge  improvement.  Several  images  are  shown  in 
the  following  hgures  to  demonstrate  the  visual  improvement  realized  by  application 
of  the  outlier  detection  algorithm. 
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5. 6. 1  Resolution  Bar  Target  Data.  Typical  weighted  versus  unweighted  res¬ 
olution  bar  imagery  is  characterized  by  the  improvement  noted  in  Fig.  5.5.  Five,  300- 
image  sets  of  resolution  bar  imagery  were  processed  by  the  outlier  detection  algorithm. 
The  300  image  sets  were  partitioned  into  50  frame  ensembles  prior  to  introduction 
to  the  outlier  detection  algorithm.  Due  to  the  10  Hz  frame  collection  rate,  the  gen¬ 
eration  of  each  ensemble  required  5  seconds,  over  which  the  atmospheric  conditions 
were  assumed  statistically  stationary.  The  actual  seeing  condition  was  estimated  from 
the  unweighted  average  image.  An  effeetive  tq  was  derived  from  the  binary  weighted 
average  images,  and  may  be  considered  a  metric  for  image  improvement  between  the 
unweighted  and  weighted  average  images.  Table  5.3  documents  the  number  of  suspect 
frames  deleted  from  the  weighted  average  image,  as  well  as  the  improvement  noted  in 
spatial  resolution  using  knife-edge  OTF  estimation  techniques. 


Table  5.3:  This  table  describes  the  results  of  weighted  frame  averaging  5  distinct 
sets  of  resolution  bar  imagery  collected  on  an  experimental  optics  range.  Values  for 
number  of  frames  rejected,  actual  and  effective  tq  are  averaged  over  the  six  50  frame 
ensembles  within  each  300  image  data  set. 


Target  Set 

Frames  Rejected 

Knife-Edge  tq 

Effective  tq 

ro  Increase 

Set  1 

19.6 

3.9 

4.3 

10% 

Set  2 

19.2 

4.1 

4.6 

12% 

Set  3 

17.9 

4.6 

4.9 

7% 

Set  4 

24.5 

3.2 

3.6 

13% 

Set  5 

22.8 

3.6 

4.1 

14% 

5.6.2  Taetieal  Image  Datasets.  To  evaluate  the  utility  of  the  algorithm  on 
collections  of  image  frames  of  nominal  tactical  scenes,  several  sets  of  imagery  were 
processed.  Three  representative  datasets  are  presented  in  the  following  sections.  Spa¬ 
tial  resolution  improvement  was  indicated  by  the  resulting  images,  and  is  quantihed 
by  estimation  of  an  effective  seeing  condition  for  each  the  weighted  and  unweighted 
average  imagery.  Table  5.4  summarizes  the  improvement  realized  by  application  of 
the  outlier  detection  algorithm  to  these  diverse  datasets. 

5-28 


Unweighted  Average  image 


Weighted  Average  image 


(a)  Unweighted  M-60  average  image.  (b)  Binary  weighted  M-60  average  image. 


Figure  5.9:  Comparison  of  the  unweighted  (a)  and  binary  weighted  average  image 
(b)  for  a  M-60  armored  tank.  Contrast  is  slightly  improved,  with  marked  improvement 
in  spatial  resolution.  Blind  estimation  of  tq  increased  from  2.7  cm  for  the  unweighted 
image,  to  3.3  cm  after  image  weighting.  24  of  50  frames  were  deleted  from  the  ensemble 
prior  to  computing  the  average  image  shown  in  subimage  (b). 

5. 6. 2.1  M-60  Armored  Tank  Vehicle.  Fig  5.9  illustrates  the  compar¬ 
ison  of  weighted  and  unweighted  average  imagery  obtained  from  remote  imaging  of 
an  M-60  armored  tank  at  10  kilometers.  Spatial  resolution  is  markedly  increased  by 
deletion  of  outlier  frames,  although  image  contrast  is  not  appreciably  increased.  24 
of  50  frames  were  selected  for  elimination  by  the  algorithm,  yielding  a  signihcant  in¬ 
crease  in  effective  ro  for  the  imaging  conditions  from  2.7  cm  to  3.3  cm.  Figure  5.10 
shows  the  increase  in  effective  seeing  conditions  as  estimated  by  the  anisoplanatic 
blind  deconvolution  algorithm  described  in  Chapters  III  and  IV.  The  peak  of  the 
likelihood  vs.  ro  curve  is  signihcantly  shifted  to  the  right  by  removal  of  suspected 
outlier  frames. 

5. 6. 2. 2  M-923  5-Ton  Truck  with  Structures.  Fig  5.11  shows  a  repre¬ 
sentative  scene  composed  of  a  military  5-ton  truck  against  a  background  of  a  building 
and  a  water  tower.  Signihcant  vertical  registration  blur  is  evident  by  inspection  of 
the  vehicle  headlights.  The  algorithm  automatically  discarded  21  of  the  50  frames  in 
the  ensemble,  producing  an  average  image  that  has  reduced  blur  and  dramatically  en- 
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Maximum  Likelihood  r  is  0.027  aftef  185  iterations  Maximum  Likelihood  r  is  0.033  after  1 51  iterations 


(a)  Unweighted  M-60  tq  plot.  (b)  Binary  weighted  M-60  tq  plot. 


Figure  5.10:  Comparison  of  the  estimated  effective  seeing  conditions  between  the 

unweighted  (a)  and  binary  weighted  (b)  average  images.  The  peak  in  the  maximum 
likelihood  estimator  occurs  at  tq  =  2.7  cm  for  the  unweighted  average,  and  increases 
to  ro  =  3.3  for  the  case  of  the  binary  weighted  average  image. 

hanced  contrast.  Blind  deconvolution  estimation  of  the  seeing  condition  has  increased 
from  approximately  4.2  to  4.8  cm,  indicating  notable  improvement  the  spatial  resolu¬ 
tion  of  the  hnal  image. 

5. 6. 2. 3  Scud  Missile  Imagery.  Fig  5.12  shows  a  Scud  missile  transport- 
erector-launch  (TEL)  vehicle  against  a  fairly  nondescript  background.  Application  of 
the  outlier  removal  technique  improved  the  image  contrast,  and  slightly  but  notably 
reduced  the  image  blur.  The  effective  tq  was  increased  by  more  than  11%  from 
3.5  cm  to  3.9  cm.  Slightly  better  wheel  dehnition  and  edge  sharpness  is  apparent. 
Inspection  of  the  individual  frames  of  this  dataset  revealed  that  several  frames  were 
poorly  registered  by  the  fast  vector  correlation  algorithm,  and  that  some  of  the  frames 
suffered  from  slight  anisoplanatic  warping  effects. 

5.7  Conclusions  and  Discussion 

Accurate  identihcation  of  outlier  image  frames  aids  the  image  processor  by  re¬ 
ducing  image  blur  due  to  averaging  frames  of  a  collected  image  ensemble.  The  max- 
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Figure  5.11:  Comparison  of  the  unweighted  (a)  and  binary  weighted  average  image 
(b)  for  a  M-923  5-ton  truck  and  surrounding  structures.  Note  the  increased  image 
contrast  indicated  by  the  maximum  photon  count  on  the  intensity  bars,  as  well  as 
reduced  blur  noted  by  inspection  of  the  specular  returns  from  the  vehicle  headlights. 
Blind  estimation  of  ro  increased  from  4.2  cm  for  the  unweighted  image,  to  4.8  cm  after 
image  weighting.  21  of  50  frames  were  deleted  from  the  ensemble  prior  to  computing 
the  average  image  shown  in  subimage  (b). 


3000 

2800 

2600 

2400 

2200 

2000 

1800 

(a)  Unweighted  Scud  average  image.  (b)  Binary  weighted  Scud  average  image. 


Figure  5.12:  Comparison  of  the  unweighted  (a)  and  binary  weighted  average  im¬ 

age  (b)  for  a  Scud  transporter-erector-launch  (TEL)  vehicle  against  a  non-descript 
background.  Note  the  slightly  increased  image  contrast,  as  well  as  reduced  blur  noted 
by  inspection  of  the  horizontal  edges  of  the  TEL  and  wheels.  Blind  estimation  of  tq 
increased  from  3.5  cm  for  the  unweighted  image,  to  3.9  cm  after  image  weighting.  22 
of  50  frames  were  deleted  from  the  ensemble  prior  to  computing  the  average  image 
shown  in  subimage  (b). 
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Table  5.4:  Table  describing  the  results  of  weighted  frame  averaging  3  diverse  sets  of 
tactical  target  imagery  collected  on  an  experimental  optics  range.  Values  for  number 
of  frames  rejected,  actual  and  effective  tq  are  calculated  for  each  of  the  three  50  frame 
ensembles. 


Target  Set 

Frames  Rejected 

Actual  ro 

Effective  ro 

ro  Increase 

M-60  Tank 

24 

2.7 

3.3 

20% 

M-923  Truck 

21 

4.2 

4.8 

14% 

SCUD  TEL 

22 

3.5 

3.9 

11% 

imum  likelihood  framework  derived  in  Sec.  5.3  appears  to  be  a  highly  effective  tool 
for  identihcation  of  such  frames. 

For  the  case  of  simulated  imagery,  the  technique  appears  to  be  highly  effective 
for  cases  of  relatively  poor  seeing  condition,  and  less  effective  as  the  seeing  parameter 
approaches  the  optical  aperture  diameter.  This  result  is  satisfying,  as  the  formation  of 
anisoplanatically  warped  images  is  less  likely  for  better  seeing  conditions,  and  global 
tip  and  tilt  are  dramatically  reduced.  In  cases  of  poor  seeing  conditions,  the  rela¬ 
tive  frequency  of  capturing  a  frame  such  as  the  warped  image  noted  in  Fig.  5.7(b) 
becomes  more  probable.  For  imagery  collected  during  seeing  conditions  character¬ 
ized  by  vq  sizes  approaching  the  optical  aperture,  the  technique  appears  unnecessary. 
There  remains  the  possibility  of  mis-registration  of  the  imagery  due  to  false  correla¬ 
tion  peaks  specihc  to  some  datasets.  This  effected  was  noted  during  the  processing 
of  experimentally  collected  data,  but  was  not  encountered  during  reduction  of  the 
simulated  data,  perhaps  due  to  the  simplicity  of  the  resolution  pattern  used  in  the 
study. 

The  results  obtained  from  processing  experimentally  collected  data  were  very 
promising.  The  concept  of  effective  tq  was  introduced  in  Sec.  5.6.1  to  demonstrate  the 
quantitative  improvement  of  the  average  images  formed  by  binary  frame  weighting 
of  the  component  ensemble  imagery.  The  increase  in  effective  seeing  conditions  was 
notable  in  all  imagery  processed.  The  effect  was  dramatic  for  the  cases  where  false 
correlation  peaks  prevented  accurate  registration  of  images.  As  an  example,  the 
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imagery  used  to  form  the  average  image  shown  in  Fig.  5.5  was  populated  by  several 
frames  that  did  not  register  accurately  due  to  a  bright  specular  return  from  the  chrome 
door  handle  of  the  vehicle  that  supported  the  resolution  target  board.  This  return 
was  erroneously  correlated  with  a  portion  of  the  resolution  target  board  in  some  image 
frames.  The  detection  of  these  frames  is  clear  and  distinct  as  demonstrated  by  the 
plot  of  likelihood  ratios  shown  in  Fig.  5.13. 


Likehood  Ratio  as  a  Function  of  Suspect  Frame  Number 


Figure  5.13:  Likelihood  ratios  calculated  for  a  set  of  experimentally  collected  im¬ 
agery  of  a  supported  resolution  target  board.  The  5  notable  outlier  frames  correspond 
to  images  that  were  not  accurately  registered  to  the  ensemble  mean  due  to  false  cor¬ 
relation  peaks  from  a  bright  specular  return. 

The  imagery  of  the  M-60  tank  discussed  in  Sec.  5.6.2. 1  is  an  example  of  data 
corrupted  by  both  anisoplanatic  warping  and  the  general  inability  of  the  registration 
algorithm  to  accurately  align  the  image  frames  given  the  severe  speckle  noise  of  the 
data.  Distinct  features  begin  to  become  apparent  in  the  resulting  image  after  outlier 
removal,  including  the  enhancement  of  the  gun  barrel  and  dehnition  of  the  tread 
wheels.  Such  details  become  important  to  the  tactical  operator  when  faced  with  the 
tasks  of  vehicle  identihcation  or  distinction  between  friend  and  foe.  The  plots  shown 
in  Fig.  5.10  depict  a  typical  output  from  the  blind  deconvolution  algorithm  described 
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in  detail  in  Chapters  III  and  IV.  The  increase  in  the  effective  seeing  condition  of  more 
than  20%  is  signihcant.  The  images  of  Figures  5.11  and  5.12  show  similar  contrast 
and  sharpness  enhancement,  and  demonstrate  effective  seeing  condition  increases  of 
14%  and  11%  respectively. 

5. 8  Summary 

A  novel  Bayesian  technique  was  developed  to  identify  those  frames  within  an 
ensemble  of  coherently  collected  imagery  that  tended  to  reduce  the  likelihood  of  the 
composite  ensemble.  It  appears  that  the  choice  of  the  maximum  likelihood  cost  func¬ 
tion  is  useful  under  a  framework  where  frames  corresponding  to  low  likelihoods  reduce 
the  contrast  and  spatial  frequency  content  of  the  ensemble  averaged  image. 

An  initial  research  effort  attempted  to  iteratively  solve  for  the  individual  con¬ 
tinuous  weights  of  each  image  that  maximized  the  overall  likelihood  of  the  entire 
ensemble.  Under  the  assumption  of  individual  pixel  distributions,  as  well  as  statisti¬ 
cal  independence  of  the  collected  images,  a  composite  likelihood  of  the  image  ensemble 
was  derived  using  the  negative  binomial  statistical  model  for  detected  intensity  of  a 
coherent  imaging  system.  It  was  found  that  the  iterative  algorithm  suffered  from 
numerical  expense.  Additionally,  it  was  difficult  to  determine  a  suitable  termination 
criteria.  Due  to  these  limitations,  a  simpler  model  was  developed  to  investigate  the 
applicability  of  a  binary  frame  weighting  model,  whereby  frames  assigned  a  weight  of 
unity  were  retained,  while  those  assigned  a  weight  of  zero  were  discarded.  A  likelihood 
ratio  test  provided  a  convenient  and  expeditious  mechanism  to  assign  these  binary 
weights. 

Application  of  the  binary  weights  to  form  a  weighted  frame  average  resulted  in 
ensemble  average  images  that  had  improved  contrast  and  spatial  resolution,  as  fur¬ 
ther  indicated  by  improvement  of  the  apparent  seeing  conditions  through  which  the 
individual  images  were  formed.  The  applications  of  this  system  are  numerous,  and 
are  not  limited  to  coherent  imaging  systems  as  a  slight  modihcation  of  the  negative 
binomial  distribution  to  the  Poisson  distribution  yields  an  estimator  for  incoherent 


5-34 


image  outlier  detection.  Any  imaging  system  faced  with  the  problems  of  anisopla- 
natism  and  tip-tilt  removal  would  beneht  from  the  selective  removal  of  images  that 
do  not  positively  contribute  to  the  frame  average.  Such  applications  include  tactical 
airborne  or  ground-based  wide  FOV  imaging,  as  well  as  astronomical  ground-based 
imaging  systems  that  process  large  ensembles  of  short-exposure  imagery. 
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VI.  Conclusions  and  Summary 


Restoration  of  images  collected  from  the  backscattered  emissions  of  a  remote 
target  illuminated  by  partially  coherent  illumination  is  a  challenging  problem 
that  offers  system  designers  signihcant  benehts  over  traditional  imaging  systems.  Al¬ 
though  passive  infrared  (IR)  imaging  systems  operating  in  the  thermal  region  offer 
the  user  moderate  resolution  without  the  need  for  active  illumination  of  the  scene, 
reliance  on  ambient  illumination  often  presents  operational  difficulties.  For  passive 
IR  imagers,  the  ambient  thermal  contrast  ratio  presents  difficult  challenges  for  ap¬ 
plication  to  target  detection  and  recognition  during  the  periods  of  thermal  crossover 
that  occur  during  dusk  and  dawn.  In  addition,  such  systems  rely  on  emissions  in  the 
8-12  micron  wavelength  region.  Such  wavelengths  are  an  order  of  magnitude  longer 
than  typical  high-power  laser  illuminators  based  on  NdiYag  technology  which  operate 
in  the  1  to  1.5  micron  region.  Since  the  resolving  power  of  an  optical  imaging  system 
follows  an  inverse  linear  relation  with  wavelength,  roughly  an  order  of  magnitude  of 
resolution  is  gained  by  reconstructing  imagery  at  the  shorter  wavelengths  typical  of 
high-power  solid-state  laser  illuminators.  In  addition  to  the  enhancement  of  basic 
resolution  limits,  the  use  of  active  illumination  introduces  tactical  and  strategic  flex¬ 
ibility  impossible  with  passive  incoherent  illumination.  Gated  laser  vision  systems 
allow  reduction  of  noise  by  way  of  accurate  shutter  control  in  unison  with  beam  pulse 
timing.  Such  a  system  allows  tremendous  flexibility  in  the  elimination  of  visible  clut¬ 
ter.  For  example,  targets  obscured  by  camouflage  netting  may  be  better  resolved 
by  first  gating  out  the  camouflage  noise  and  operating  on  only  the  data  within  the 
tightly  gated  region  surrounding  the  target  depth.  In  addition,  such  mechanization 
provides  accurate  range  information  to  different  regions  of  the  target  held,  although 
the  challenges  associated  with  this  application  were  not  investigated  in  this  research 
effort. 
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6.1  Summary  of  MAP  Estimation  of  Partially  Coherent,  Anisoplanat- 
ically  Distorted  Imagery 

The  accurate  restoration  of  images  created  by  an  active  coherent  vision  system 
comes  at  the  expense  of  some  particular  difficulties.  Even  without  the  deleterious 
effects  of  atmospheric  turbulence,  the  coherent  nature  of  the  active  illumination  causes 
highly  speckled  imagery  that  is  often  unsuitable  for  presentation  to  the  operator 
without  some  form  of  image  post-processing.  The  research  described  in  the  preceding 
chapters  presupposes  that  individual  images  are  first  averaged  to  reduce  the  gross 
effects  of  laser  as  well  as  atmospheric  speckle.  Whether  the  image  processor  uses  a 
single,  or  an  ensemble  average  of  partially  coherent  imagery,  the  accurate  formulation 
of  a  likelihood-based  image  estimator  depends  hrmly  on  the  underlying  assumption 
of  the  probabilistic  distribution  of  the  detected  illumination.  The  negative  binomial 
probability  mass  function  has  been  demonstrated  to  be  a  very  accurate  model  that 
conveniently  extends  from  fully  developed  speckle  imagery,  to  images  formed  from 
relatively  incoherent  lasers  or  even  incoherent  illumination  in  the  more  extreme  case. 
In  this  regard,  the  research  described  here  is  easily  extended  along  the  continuum 
of  coherency  ranging  from  laser  illuminators  with  extremely  long  coherence  times 
to  passively  illuminated  scenes.  The  latter  case  is  merely  a  convenience  due  to  the 
extension  of  the  negative  binomial  distribution  to  the  Poisson  distribution  for  the 
limiting  situation  of  very  large  speckle  parameters. 

Atmospheric  turbulence  causes  tremendous  image  distortion  for  scenarios  where 
long  slant-range  paths  over  low-altitude  turbulence  is  unavoidable.  Turbulence  close 
to  the  optical  aperture  causes  the  majority  of  phase  abberation  in  the  detected  im¬ 
agery.  Unlike  satellite  space- vehicle  imaging  systems,  airborne  or  ground-based  imag¬ 
ing  systems  will  suffer  high  levels  of  atmospheric  distortion.  Adaptive  optic  wave- 
front  pre-distortion  allows  highly  effective  image  restoration  for  cases  where  the  at¬ 
mospheric  distortion  may  be  estimated  in  near  real-time.  However,  such  systems  are 
large,  expensive  and  computationally  and  mechanically  complex.  Application  of  AO 
technology  to  space  and  power-constrained  platforms  is  a  rich  area  of  research  that 
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will  bear  fruit  in  the  coming  decades.  In  the  meantime,  there  is  a  growing  need  for 
post-processing  algorithms  that  can  effectively  mitigate  the  effects  of  turbulence  on 
images  formed  over  long  distances  through  large  volumes  of  the  atmosphere.  To  fur¬ 
ther  complicate  the  problem,  the  FOV  of  candidate  tactical  systems  is  relatively  wide, 
especially  compared  to  those  of  the  astronomical  research  community,  where  much  of 
the  research  in  image  restoration  algorithms  has  been  conducted  over  the  past  sev¬ 
eral  decades.  Effective  image  restoration  algorithms  must  deal  with  the  high  levels  of 
anisoplanatism  that  occur  in  these  systems,  even  in  conditions  where  the  turbulence 
is  relatively  moderate. 

Chapters  III  and  IV  built  upon  the  frameworks  presented  in  Chap.  II  to  build 
maximum  a  posteriori  estimators  for  cases  where  the  imaging  system  was  considered 
spatially-invariant  and  spatially- variant  respectively.  Several  sources  of  blur  conspire 
to  reduce  high  spatial  frequency  detail  in  the  hnal  detected  image.  For  moderately 
turbulent  conditions  over  long  horizontal  or  slant  paths,  the  atmospheric  seeing  condi¬ 
tion,  parameterized  by  Fried’s  seeing  parameter  tq,  becomes  much  more  of  a  limiting 
system  factor  than  the  limits  imposed  by  the  physical  aperture.  Tip  and  tilt  com¬ 
ponents  of  the  atmospheric  random  phase  delays  cause  signihcant  blur  due  to  linear 
translation  in  the  orthogonal  axes  of  the  image.  Fortunately,  by  the  judicious  use  of 
robust  registration  algorithms,  most  if  not  all  of  this  motion  blur  may  be  effectively 
removed  from  the  resulting  ensemble  average  image.  A  more  difficult  problem  is  en¬ 
countered  when  attempting  to  remove  the  blur  caused  by  the  uncorrelated  motion 
of  the  many  isoplanatic  sub-image  patches  that  occur  due  to  the  shift-variant  nature 
of  wide  FOV  optical  systems  operating  through  even  moderate  levels  of  atmospheric 
turbulence.  AO  systems  with  multiple  points  of  reference  or  “guide  stars”  are  effective 
tools  to  estimate  and  deconvolve  the  spatially  variant  OTF  that  describe  this  process. 
The  simpler  approach  described  in  Chap.  IV  provides  a  means  to  estimate  the  addi¬ 
tional  average  blur  created  in  the  average  image  due  to  the  spatially  variant  imaging 
system.  Given  such  a  parameterized  model  for  anisoplanatic  blur,  the  deconvolution 
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kernel  of  the  MAP  estimator  presented  in  Chap.  Ill  is  improved  to  better  process 
imagery  collected  through  a  wide  FOV  system. 

Chapter  V  explored  the  additional  benefit  gained  in  effective  seeing  condition 
enhancement  by  the  selective  removal  of  ensemble  image  frames  from  the  averaged 
image  introduced  to  the  deconvolution  algorithm.  Although  such  processing  intro¬ 
duces  a  bias  in  the  estimation  of  the  actual  seeing  condition  parameter,  the  goal  of 
obtaining  higher  spatial  frequency  content  within  the  deconvolved  image  was  attained. 
Identihcation  and  removal  of  suspect  ensemble  images  is  a  convenient  way  to  enhance 
imagery  collected  from  a  system  with  relatively  high  frame  sample  periods.  As  an 
example,  a  system  with  a  moderately  fast  frame  rate  might  discard  25  of  50  frames  in 
an  ensemble.  If  50  frames  are  deemed  necessary  to  satisfy  image  SNR  requirements 
due  to  low  photon  counts  resulting  from  long  range  imaging  scenarios,  an  additional 
50  frames  might  be  collected  during  some  acceptably  short  dwell  period. 

6.2  Research  Contribution  Summary 

Several  specihc  and  signihcant  research  contributions  result  from  the  work  dis¬ 
cussed  in  this  document.  These  contributions  are  intended  for  application  to  the  field 
of  coherent  image  restoration,  but  may  be  extended  to  the  broader  held  of  incoherent 
illumination  in  many  circumstances. 

6.2.1  Restoration  of  Remote  Seene  Imagery  Illuminated  by  Partially  Coherent 
Light.  The  accurate  restoration  of  imagery  captured  by  a  partially  coherent  laser 
vision  system  is  hampered  by  the  speckle  that  is  caused  by  the  physical  interaction 
of  the  illuminating  beam  with  the  target  surface,  as  well  as  the  speckle  created  by 
the  random  delays  imposed  by  the  turbulent  atmosphere  between  the  optical  system 
and  the  target.  The  physically  based  propagation  model  of  Chap.  II  provides  a  sta¬ 
tistical  means  by  which  the  detected  intensity  of  a  coherently  illuminated  target  may 
be  reconstructed  via  maximum  likelihood  estimation  techniques.  This  research  led 
to  the  development  of  such  an  estimator  in  Chap.  Ill  to  jointly  estimate  the  remote 
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scene  together  with  the  actual  seeing  conditions  under  which  the  scene  was  imaged, 
parameterized  by  tq.  The  utility  of  such  an  estimator  is  quite  general.  In  the  limit 
as  the  speckle  parameter  is  allowed  to  grow  large,  the  imaging  situation  closely  re¬ 
sembles  incoherent  imaging.  In  such  cases,  the  blind  deconvolution  MAP  estimation 
technique  might  allow  the  accurate  estimation  of  atmospheric  seeing  conditions  over 
long  optical  paths  without  the  use  of  expensive  scintillometry  equipment  or  special 
illuminator  sources.  A  simple  experimental  observation  of  a  scene  with  sufficiently 
high  spatial  detail  might  be  all  that  is  required  to  yield  accurate  estimates  of  the  com¬ 
posite  horizontal  or  slant-path  integrated  turbulence  between  the  target  and  optical 
system. 

6.2.2  Anisoplanatic  OTF  Describing  Wide  FOV  Imaging  Systems.  For  wide 
FOV  systems,  the  absolute  level  of  turbulence  typical  of  terrestrial  operating  scenarios 
causes  highly  anisoplanatic  viewing  conditions.  The  additional  blur  that  arises  due  to 
the  spatially  variant  OTF  must  be  properly  accounted  for.  The  research  of  Chap.  IV 
provided  a  concise  and  effective  description  of  the  quantitative  effects  of  this  type  of 
optical  degradation,  and  presented  a  compact  model  for  the  parameterization  of  these 
effects.  When  incorporated  into  the  MAP  blind  deconvolution  algorithm  described 
initially  in  Chap.  Ill,  the  model  was  found  to  better  represent  the  true  atmospheric 
conditions  used  to  image  the  remote  scene.  The  extension  to  anisoplanatic  seeing 
conditions  enhances  the  utility  of  a  calibrated  seeing  condition  estimator.  For  tur¬ 
bulence  levels  that  cause  wide  FOV  optical  systems  to  exceed  the  isoplanatic  angle, 
the  monitor  remains  effective  in  presenting  seeing  condition  estimates  not  appreciably 
biased  by  the  introduction  of  spatially  variant  blur  into  the  ensemble  images  collected 
by  the  system. 

6.2.3  Seeing  Condition  Monitor.  Although  much  of  the  motivation  for  the 
development  of  the  algorithms  presented  in  this  research  lies  in  the  ability  of  the  im¬ 
age  restoration  process  to  enhance  spatial  resolution  and  image  detail,  an  important 
side-effect  is  the  accurate  estimation  of  the  atmospheric  conditions  under  which  the 
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images  were  collected.  A  distinct  contribution  is  made  to  those  researchers  who  seek 
new  methods  to  accurately  estimate  seeing  conditions  through  a  variety  of  different 
and  evolving  atmospheres.  As  an  example,  a  compact  seeing  monitor  may  be  con¬ 
structed  whereby  a  remote  scene  is  illuminated  by  partially  coherent  light,  and  the 
techniques  and  algorithms  described  in  Chapters  III  and  IV  are  used  to  provide  ac¬ 
curate  estimates  for  Fried’s  seeing  parameter.  The  new  method  is  valuable  due  to 
non-reliance  on  particular  target  scenery,  as  might  be  required  by  competing  systems. 
As  will  be  discussed  in  Sec.  6.3.3,  the  estimation  algorithms  may  be  extended  to  pro¬ 
cess  incoherently  illuminated  scenes,  such  as  those  illuminated  by  ambient  light.  Such 
extension  is  due  to  the  generality  of  the  model  used  to  quantify  the  statistics  of  the 
detected  light  of  the  imaging  device. 

6.2.4  Outlier  Detection  and  Binary  Weighted  Frame  Averaging  of  Ensembles 
of  Coherently  Detected  Imagery.  Chapter  V  introduced  the  utility  of  an  algorithm 
that  seeks  improved  spatial  resolution  by  the  automatic  selection  of  particular  frames 
within  an  ensemble  of  images  collected  through  a  random  atmosphere.  The  sources 
of  image  corruption  included  the  basic  speckle  mechanism  described  above,  but  also 
the  random  turbulence  and  mis-registration  between  images  within  the  ensemble.  Of 
course,  such  a  system  would  tend  to  produce  optimistically  biased  estimates  of  the 
seeing  condition  when  employed  as  a  seeing  condition  monitor,  however,  the  system  is 
of  great  value  to  communities  that  wish  to  regain  high  spatial  frequency  information 
from  a  series  of  averaged  images.  For  example,  given  a  laser  vision  system  with 
sufficiently  high  frame  rate,  large  numbers  of  speckle  images  may  be  collected  over 
acceptably  brief  dwell  periods.  Under  these  circumstances,  it  may  be  more  efficient 
to  simply  discard  frames  that  are  heavily  corrupted  by  turbulence,  or  difficult  to 
accurately  register  due  to  anisoplanatic  warping,  rather  than  spend  inordinately  large 
amounts  of  computing  resources  in  an  attempt  to  repair  the  corrupted  images  by 
using  image  de-warping  techniques  or  alternate  registration  algorithms. 


6-6 


6.2.5  Speckle  Parameter  Estimator.  The  application  of  the  MAP  estima¬ 
tion  techniques  discussed  above  rely  on  fairly  accurate  estimation  of  the  laser  speckle 
parameter  Ai  that  characterizes  the  NB  distribution  of  the  detected  image  intensity. 
The  speckle  parameter  is  not  hxed,  and  depends  on  complex  relationships  between 
the  laser  coherence  time,  the  detector  gating  period,  the  amount  of  laser  beam  scin¬ 
tillation,  and  perhaps  other  variables  that  are  difficult  to  measure  directly.  Accurate 
estimation  of  Ai  as  described  in  Sec.  2.5  provides  an  information  theoretic  approach 
to  the  calculation  of  the  effective  speckle  parameter  that  parameterizes  the  random 
process  assumed  to  govern  the  MAP  estimators  of  Chapters  III  and  IV. 

6.2.6  Effects  of  Image  Quantization  and  Scaling.  The  effects  of  image  inten¬ 
sity  quantization  and  scaling  were  recognized  in  the  early  stages  of  this  research  and 
caused  some  difficulty  in  the  interpretation  of  the  simulated  and  measured  results. 
Section  2.7  discussed  the  investigation  and  effects  produced  by  these  phenomena. 
Fielded  image  processing  systems  will  likely  involve  compromises  involving  image 
quantization.  In  most  cases,  accurate  image  scaling  must  be  performed  to  relate  the 
recorded  image  intensities  to  the  number  of  photons  received.  This  relationship  is 
important  because  the  MAP  estimators  were  constructed  using  statistical  models  of 
photon  arrival  at  the  detector.  Without  attention  to  image  scaling  and  quantization, 
the  measured  data  will  not  follow  the  expected  statistics  of  the  imaging  system.  This 
research  presents  the  tools  necessary  by  which  the  effects  of  image  quantization  and 
scaling  may  be  understood  in  the  context  of  coherent  image  restoration  using  MAP 
estimation,  and  raw  image  data  may  be  accurately  calibrated  for  introduction  into 
these  and  similar  algorithms. 

6.3  Future  Research  Considerations 

Although  the  research  described  in  this  document  represents  a  fairly  complete 
treatment  of  image  restoration  of  partially  coherent  remote  imagery,  several  research 
efforts  might  extend  the  utility  of  the  described  methods. 
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6. 3. 1  Speed  Improvement  of  Blind  Deeonvolution  Algorithm  for  Real-Time  Ap- 
plieations.  Although  the  blind  deconvolution  algorithm  described  in  Chap.  Ill  con¬ 
verges  fairly  rapidly  to  approach  the  maximum  likelihood  estimate  of  the  scene  for  a 
particular  value  of  tq,  the  space  of  values  for  tq  is  fairly  large,  requiring  a  search  over 
the  entire  space.  Due  to  the  apparent  monotonicity  of  the  likelihood  curve  for  each 
value  of  ro,  it  may  prove  feasible  to  obtain  a  coarse  estimate  of  the  seeing  condition 
using  only  several  iterations  at  each  value  in  rg-space.  The  likelihood  of  the  scene 
may  then  be  maximized  only  for  neighboring  values  of  tq. 

Although  this  approach  appears  promising,  a  much  faster  algorithm  might  be 
realized  if  all  calculations  could  be  performed  in  the  frequency  domain,  where  convolu¬ 
tions  are  simply  evaluated  using  circular  convolution  by  way  of  fast-Fourier  transform 
techniques.  Referring  to  Eqn.  3.15,  repeated  here  for  convenience. 
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it  is  clear  that  both  the  numerator  and  the  denominator  of  the  right-hand  side  may 
be  evaluated  using  Fourier  domain  convolution.  However,  at  each  iteration  step, 
the  right-hand  side  fraction  must  be  transformed  back  to  the  spatial  domain  in  or¬ 
der  to  pointwise  multiply  by  the  previous  scene  estimate.  This  process  requires  the 
2-dimensional  Fourier  transformation  of  a  fairly  large  matrix.  A  more  efficient  im¬ 
plementation  would  be  realized  if  the  entire  operation  could  be  cast  into  the  Fourier 
domain,  with  the  requirement  to  return  to  the  spatial  domain  only  twice,  at  the  ini¬ 
tialization  and  end  of  the  iterative  process.  However,  the  negative  binomial  statistics 
of  the  detected  intensity  presents  a  difficult  analytical  problem.  A  tempting  approach 
would  be  to  represent  the  negative  binomial  statistics  into  a  more  easily  transformable 
distribution,  such  as  the  Gaussian  PDF. 

The  study  of  Gaussian  mixture  models  (GMM)  provides  a  highly  useful  tool 
that  has  been  applied  to  many  statistical  modeling  applications  [5] .  GMMs  allow  the 
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fairly  arbitrary  representation  of  virtually  any  statistical  distribution  with  one  or  more 
Gaussian  basis  functions.  Since  the  Fourier  operator  is  linear,  the  negative  binomial 
statistics  of  the  coherent  photodetection  process  might  be  well  approximated  by  some 
number  of  Gaussian  basis  functions,  and  easily  recast  into  the  Fourier  domain  to  allow 
complete  reconstruction  of  the  image  without  alternating  between  both  domains. 

This  course  of  research  might  lead  to  an  extremely  fast,  reliable  and  robust 
algorithm  with  signihcant  practical  value  to  researchers  and  system  designers  requiring 
fast,  accurate  deconvolution  of  a  wide  range  of  coherently  imaged,  turbulence  degraded 
wide  FOV  scene  data. 

6.3.2  Proof  of  Convergence  of  the  Iterative  Algorithms.  A  further  beneht 
of  the  application  of  GMMs  to  this  problem  might  yield  a  formulation  of  the  decon¬ 
volution  algorithm  in  light  of  the  optimality  of  the  expectation-maximization  (EM) 
algorithm  [5,50].  Although  the  coherent  blind  deconvolution  algorithm  is  derived  in 
a  maximum  a  posteriori  framework,  it  has  not  been  proven  that  the  algorithm  in¬ 
creases  the  likelihood  as  iterations  progress,  despite  all  indications  that  support  this 
conclusion.  One  of  the  exciting  properties  of  the  application  of  the  EM  algorithm  is 
that  “it  has  proved  to  be  a  valuable  tool  for  many  problems,  since  it  provides  an  ele¬ 
gant  approach  to  bypass  difficult  optimization  and  integrations  required  in  Bayesian 
estimation  problems”  [50].  The  difficulty  of  applying  the  EM  algorithm  frequently 
occurs  during  the  formulation  of  the  E-step,  where  the  conditional  densities  of  the 
hidden  variables  must  be  determined.  The  transformation  of  statistics  from  negative 
binomial  to  Gaussian  using  GMMs  may  make  such  an  approach  tractable. 

6.3.3  Extension  to  Incoherently  Collected  Imagery.  The  techniques  de¬ 
scribed  in  the  research  may  be  extended  from  partially  coherent  to  relatively  inco¬ 
herent  illumination,  as  is  encountered  in  the  vast  majority  of  imaging  scenarios.  The 
extension  is  a  natural  result  of  the  degree-of-freedom  introduced  by  the  speckle  pa¬ 
rameter  M.  in  the  negative  binomial  distribution  that  describes  the  photon  count  of 
partially  coherent  illumination  at  the  imaging  detector.  As  an  example  of  a  diverse 
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application  of  this  research,  smaller  astronomical  observatories  might  restore  images 
of  remote  objects  that  are  subject  to  the  same  ill  effects  of  terrestrial  objects  when 
viewed  through  large  aperture  systems  limited  by  the  atmospheric  coherence  diam¬ 
eter.  A  common  example  is  the  observation  of  the  Earth’s  Moon  or  the  surface  of 
the  Sun.  Both  objects  have  relatively  large  angular  extent  and  often  require  image 
restoration  techniques  that  effectively  cope  with  anisoplanatic  imaging  conditions. 
Automatic  frame  selection  of  such  images  provides  further  spatial  detail  that  would 
otherwise  be  difficult  to  cull  from  a  large  set  of  speckle  imagery.  Although  larger, 
better  equipped  observatories  might  have  AO  image  enhancement  capability,  such 
systems  often  require  a  portion  of  the  illumination  to  be  used  for  wavefront  estima¬ 
tion,  thus  reducing  the  final  photon  count  at  the  imaging  detector.  In  addition,  the 
added  expense  and  complexity  of  a  full  AO  system  might  not  be  justified  for  smaller 
observatory  missions. 

6.3.4  Fusion  of  Imaging  Correlography  Information  with  Imaged  Data. 
Imaging  correlography  is  an  interesting  field  of  research  that  seeks  to  reconstruct 
an  image  from  the  Fourier  modulus  of  the  fields  collected  at  the  aperture  plane,  with¬ 
out  the  requirement  of  focusing  the  field  on  an  imaging  detector  array.  Fienup  and 
others  have  reported  good  results  in  the  synthesis  of  images  obtained  by  the  coherent 
illumination  of  reflected  laser-speckle  intensitity  patterns  [36,74]. 

As  presented,  the  recovery  of  high-resolution  images  from  the  Fourier  modulus 
collected  at  the  aperture  is  a  computationally  intense  process  that  requires  many 
data  frames  to  achieve  suitably  high  SNR.  However,  the  collection  of  only  Fourier 
modulus  data  at  the  aperture  discards  phase-dominated  atmospheric  distortion  as 
evidenced  by  the  model  presented  in  Sec.  2.3.  A  potential  research  avenue  would  be 
the  exploration  of  the  advantages  of  fusing  these  data  with  the  imaged  data  collected 
by  the  detector  behind  the  aperture  lens.  Since  the  Fourier  modulus  data  are  relatively 
unaffected  by  the  phase-dominated  turbulence  variations,  it  is  expected  that  high 
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frequency  information  available  in  the  correlography  data  would  complement  the  low- 
pass  hltered  data  that  results  at  the  image  plane  due  to  the  atmospheric  turbulence. 

The  difficulty  of  this  fusion  effort  lies  in  the  ability  to  cast  the  coherent  illumi¬ 
nation  statistics  into  the  aperture  domain.  Since  the  negative  binomial  distribution 
applies  only  to  the  intensity  detection  at  the  imaging  detector  plane,  this  distribution 
is  not  necessarily  valid  at  the  aperture  plane.  As  discussed  in  Sec.  6.3.1,  the  recast 
of  the  negative  binomial  distribution  into  the  Fourier  domain  is  not  trivial.  However, 
research  might  be  conducted  to  arrive  at  a  direct  physical  model  of  the  statistics  of  the 
Fourier  modulus  of  a  scene  illuminated  by  partially  coherent  light  backscattered  from 
a  target  scene.  Armed  with  this  statistic,  the  additional  information  would  be  added 
to  the  likelihood  equation  and  maximized  using  a  similar  iterative  blind  deconvolution 
algorithm.  It  is  expected  that  the  fusion  of  the  aperture  derived  correlography  data 
with  detected  imagery  will  improve  the  MAP  estimator  algorithm  performance. 

6.4  Final  Thoughts 

Active  illumination  of  remote  targets  using  partially  coherent  laser  illumina¬ 
tion  provides  system  designers  unprecedented  levels  of  operational  freedom.  Accurate 
restoration  of  the  detected  imagery  is  essential  to  the  success  of  such  systems.  Al¬ 
though  adaptive  optics  provide  an  attractive  methodology  by  which  imagery  may  be 
effectively  enhanced,  the  helding  of  robust,  compact  and  reliable  systems  is  still  many 
years  in  the  future.  Image  post-processing  techniques  provide  an  immediate  solution 
to  a  difficult  and  rewarding  problem.  The  techniques  described  in  this  research  are 
presented  as  stepping  stones  toward  the  goal  of  realizing  useful  and  robust  image 
reconstruction  systems  for  terrestrial  imaging  scenarios  of  interest  to  a  diverse  range 
of  system  operators. 
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Appendix  A.  Direct  Solution  of  Frame  Average  Weights 


The  iterative  solution  derived  in  Sec.  5.3.1  is  computationally  expensive  for  large 
ensembles  of  large  images.  Additionally,  stopping  criteria  for  the  iterative  al¬ 
gorithm  is  difficult  to  establish.  For  these  reasons,  a  direct  solution  to  the  likelihood 
maximization  process  is  attractive.  However,  direct  maximization  of  Eqn.  5.5  ap¬ 
pears  mathematically  intractable.  An  alternative  solution  is  offered  by  taking  several 
liberties  with  the  underlying  probability  mass  function.  The  derivation  begins  with 
several  simplifying  assumptions  about  the  detected  intensity  distribution  that  admit 
Gaussian  statistics  rather  than  the  negative  binomial  distribution  of  Eqn.  5.1.  The 
Gaussian  model  provides  a  log-likelihood  function  that  is  easily  maximized  by  the 
solution  of  a  system  of  linear  equations.  Due  to  the  construction  of  the  likelihood 
function,  it  becomes  apparent  that  the  linear  system  provides  a  least-squares  solution 
to  the  problem. 


A.l  Maximizing  the  Likelihood  of  the  Weighted  Average  Ensemble 

Recall  the  log-likelihood  equation  developed  using  the  negative  binomial  statis¬ 
tics  of  the  detected  intensity  at  the  focal  plane  array. 
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where  the  mean  intensity  i{x,y)  is  formed  by 
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Maximization  of  the  likelihood  is  taken  with  respect  to  the  individual  frame 
weights  An-,  n  =  1,2,  ...,J  where  J  is  the  total  number  of  frames  in  the  ensemble. 
The  maximization  is  difficult  due  to  the  combination  of  i{x,y)  and  Ai  in  the  loga¬ 
rithms  of  the  second  and  third  terms  of  Eqn.  A.l.  On  possible  approach  is  to  assume 
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M  small  compared  to  i{x,y)  in  order  to  simplify  the  likelihood  function.  Unfortu¬ 
nately,  such  an  approximation  is  not  well  justihed  from  a  systems  approach,  and  still 
yields  an  expression  with  a  large  summation  within  the  logarithm  due  to  the  relation¬ 
ship  of  Eqn.  A. 2.  A  more  practical  approach  may  be  undertaken  by  returning  to  the 
negative  binomial  distribution  of  Eqn.  5.1.  Under  conditions  where  the  speckle  param¬ 
eter  is  fairly  high,  for  example,  Af  =  50  or  more,  the  negative  binomial  distribution 
approaches  that  of  a  Poisson  distribution,  as  illustrated  in  Fig.  2.6.  Under  moder¬ 
ately  high  photon  count  conditions,  the  Poisson  distribution  is  well  approximated  by 
a  Gaussian  distribution,  although  the  latter  admits  the  non-physical  possibility  of 
negative  intensity  values  under  low  photon  conditions.  Using  such  an  approximation, 
the  probability  of  a  pixel  of  the  image  data  given  a  pixel  of  the  average  intensity  may 
be  expressed  as 
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where  a  is  the  unspecihed  standard  deviation  of  the  noise  process,  and  the  mean  of  the 
Gaussian  distribution  is  simply  the  mean  intensity  formed  by  the  weighted  average 
according  to  Eqn.  A. 2  with  undetermined  weight  vector  A. 

For  a  particular  image  in  the  J  frame  ensemble,  the  probability  of  a  detected 
image  given  the  weighted  frame  average  is 
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and  by  assumption  of  independence  between  image  frames  within  the  ensemble,  the 
probability  distribution  for  the  entire  ensemble  becomes 
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where  d  represents  the  ensemble  of  detected  images. 
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Maximization  of  the  log-likelihood  is  mathematically  convenient,  and  Eqn.  A. 5 


may  be  expressed  in  logarithmic  format  as 
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Since  each  maximization  of  Eqn.  A. 6  must  be  calculated  with  respect  to  an 
arbitrary  frame  weight  the  hrst  term  within  the  summation  may  be  disregarded 
and  a  new  log-likelihood  may  be  written  as 
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To  maximize  this  expression,  the  derivative  with  respect  to  Aj^  may  be  calcu¬ 
lated  and  set  to  zero.  As  derived  in  Sec.  5.3.1,  the  derivative  of  the  weighted  intensity 
with  respect  to  an  arbitrary  weight  in  the  ensemble  can  be  expressed  as 
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which  allows  differentiation  of  Eqn.  A.7  to  yield 
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By  setting  the  derivative  to  zero  and  rearranging  the  order  of  summation. 
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and  since  dj^  —  djo,  y  —  is  not  a  function  of  j  in  the  ensemble  summation,  the 
relationship  may  be  rewritten  as 
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Although  Eqn.  A.  11  appears  to  be  an  unlikely  candidate  for  direct  solution, 
expressing  the  relationship  as  a  system  of  linear  equations  for  each  particular  image 
frame  in  the  ensemble  admits  solution  by  way  of  linear  algebra  techniques.  The 
left-hand-side  of  Eqn.  A.  11  is  constant  given  a  particular  choice  for  the  frame  under 
consideration,  dj^.  Let  this  constant  be  Cj^, 
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The  expression  for  the  weighted  average  intensity  of  Eqn.  A. 2  may  be  substituted 
into  the  right-hand-side  of  Eqn.  A.  11  to  yield 
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Let  the  double  summation  over  x  and  y  of  Eqn.  A.  14  represent  a  vector  of 
coefficients  Kj^  for  each  image  frame  in  the  ensemble  dj.  Under  this  framework,  the 
linear  equation  for  a  particular  frame  dj^  can  be  written  as 
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The  entire  set  of  coefficients  for  all  frames  in  the  the  ensemble  may  be  compactly 
described  by  a  square  Jx  J  matrix  K,  where  each  row  holds  a  vector  of  such  coefficients 
for  each  selected  image  in  the  ensemble.  Thus 


N  N 

Ka,b  =  -  tta,  v  -  (3^  db  -  Oib.y  -  Pb^  , 

x=l  y=l 


(A.16) 


where  a  indexes  the  row  position,  and  b  indexes  the  column  position  of  the  matrix  K. 
In  effect,  K  is  analogous  to  a  correlation  matrix  formed  by  the  pointwise  multiplication 
of  each  frame  in  the  ensemble  with  every  other  frame  in  the  ensemble.  Unlike  a 
correlation  matrix,  however,  K  is  not  normalized. 

Under  this  representation,  Eqn.  A.  11  may  be  expressed  using  vector  notation 
as  C*  =  KA  where  A  and  C  are  J  element  vectors.  If  K  is  invertible,  then  the  direct 
solution  to  the  frame  weights  may  be  found  using  the  relationship  A  =  K'^U.  The 
random  variation  of  the  images  in  the  ensemble  permit  the  inversion  of  K.  With  little 
or  no  variation  between  each  image  in  the  ensemble,  the  rank  of  K  would  clearly  be 
less  than  J .  Low  variance  between  imagery  will  increase  the  matrix  condition  and 
cause  difficulty  in  the  inversion  process.  However,  for  realistic  imaging  scenarios,  the 
elevated  image  variance  will  increase  the  rank  of  K  to  J .  Furthermore,  the  element  of 
K  are  always  positive  due  to  the  physical  detection  of  the  photon  intensity,  and  the 
matrix  is  real  symmetric.  Under  these  conditions,  a  unique  solution  for  the  weights 
may  be  determined  by  calculation  of  A  =  using  Equations.  A.  12  and  A.16. 


A. 2  Implementation  of  the  Direct  Solution 

The  direct  solution  provided  by  the  analysis  of  Sec.  A.l  minimizes  the  square 
error  between  the  weighted  average  intensity  and  the  ensemble  imagery.  It  is  not  clear 
that  such  an  estimator  for  the  weights  of  the  frame  imagery  should  improve  the  resolu¬ 
tion  of  the  estimated  average  image  by  increasing  the  high  spatial  frequency  content. 
In  practice,  it  was  found  that  despite  the  invertibility  of  the  coefficient  matrix  K, 
the  weights  calculated  using  A  =  had  only  minor  effect  on  the  weights  applied 
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to  the  ensemble  imagery.  Several  simulated  and  experimental  datasets  were  analyzed 
using  the  direct  approach,  however,  the  resulting  weighed  average  images  did  not  sub¬ 
stantially  differ  from  the  unweighted  imagery  in  terms  of  spatial  frequency  content 
or  resolution.  While  the  direct  solution  is  slightly  less  computationally  involved  in 
comparison  to  the  binary  frame  weighting  technique  discussed  in  Chap.  V,  the  lack 
of  image  improvement  suggests  that  the  application  of  a  least-squares  solution  lacks 
signihcant  merit  when  applied  to  this  particular  problem. 
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